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>The  computational  complexity  of  a  nonlinear  system  identification  technique 
is  evaluated  in  this  report.  The  identification  technique  uses  a  ‘‘black  box**'1" 
approach  requiring  measurements  only  at  system  input  and  output  terminals  and  is 
applicable  to  weakly  nonlinear  systems  whose  behavior  is  adequately  character izec 
by  a  finite  Volterra  Series. 

The  computational  aspects  of  the  technique  are  evaluated  in  terms  of  the 
complexity  of  the  calculations  and  the  complexity  of  the  system  being  implemented 
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EVALUATION 


The  process  of  system  identification  consists  of  postulating  a  valid  analytical 
mode)  for  the  system  under  consideration  and  performing  tests  on  the  system  to 
completely  specify  or  "identify"  the  parameters  which  describe  the  system 
analytical  model.  For  example,  a  linear  system  is  completely  characterized  by  its 
impulse  response,  h(t).  The  system  identification  process  for  this  linear  system 
analytical  model  consists  of  any  procedure  that  completely  determines  h(t).  The 
present  consideration  in  the  area  of  nonlinear  system  identification  is  the  deriva¬ 
tion  of  a  valid  analytical  model  for  the  nonlinear  system  under  consideration. 

The  identification  procedure  successfully  studied  is  a  black  box  technique 
where  only  input  and  output  terminal  measurements  of  the  nonlinear  system  are 
used.  The  identification  technique  is  applicable  to  a  broad  class  of  weakly 
noniinear  systems  whose  response  can  be  characterized  by  a  finite  Volterra  series. 
The  identification  procedure  involves  processing  the  input  and  output  responses  of 
a  nonlinear  system  to  obtain  a  set  of  linearly  independent  equations  which  uniquely 
define  the  parameters  of  a  functional  form  of  the  second-order  impulsa  response. 
Theoretically,  the  proposed  identification  technique  represents  a  significant  im¬ 
provement  over  existing  identification  techniques  because  of  its  black  box  formula¬ 
tion.  The  intent  of  the  study  was  to  determine  where  this  identification  technique 
can  be  practically  implemented  and  maintain  an  advantage  over  existing  tecniques. 
To  these  ends,  the  practical  implementation  constraints  have  been  developed, 
quantified  and  assessed  for  three  candidate  measurement  configurations.  The 
robustness  of  the  technique  to  nonlinear  circuits  with  many  and/or  repeated  poles 
is  the  subject  of  Part  II  of  this  final  report. 
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The  study  effort  successfully  accomplished,  in  two  parts,  a  Part  I  on 
implementation  feasibility  to  determine  the  practical  methods  and  constraints  of 
implementing  three  candidate  measurement  configurations  -  digital,  analog  and 
hybrid.  The  second  part  of  the  study  effort  successfully  focused  on  the  numerical 
computation  complexity  aspect  of  the  identification  technique  processing  to 
determine  the  class(es)  of  nonlinear  systems  for  which  the  technique  can  be 
practically  applied.  The  primary  computational  complexity  arises  from  the 
required  matrix  inversions  for  the  residue  evaluations.  Toward  the  goal  of 
alleviating  these  difficulties,  matrix  scaling,  hand  limited  approaches,  single 
exponential  inputs  (multiple  input  times)  and  dominant  pole  concepts  were  also 
developed,  quantified  and  assessed  in  the  successful  pursuit  of  the  overall  study 


Project  Engineer 


iv 


TECHNICAL  REPORT  SUMMARY  FOR  RADC-TR-79- 
NONLINEAR  SYSTEM  IDENTIFICATION  STUDY 
PART  II.  COMPUTATIONAL  COMPLEXITY  STUDY 


A.  STUDY  OBJECTIVES 

The  basic  objective  of  this  study  effort  is  to  evaluate  the 
practical  feasibility  of  a  nonlinear  system  identification  tech¬ 
nique.  The  identification  procedure  studied  is  a  black  box 
technique  where  only  input  and  output  terminal  measurements  of 
the  nonlinear  system  are  used.  The  identification  technique  is 
applicable  to  a  broad  class  of  weakly  nonlinear  systems  whose 
response  can  be  characterized  by  a  finite  Volterra  series.  The 
identification  procedure  involves  processing  the  input  and  out¬ 
put  responses  of  a  nonlinear  system  to  obtain  a  set  of  linearly 
independent  equations  that  uniquely  define  the  parameters  of  a 
functional  form  of  the  second-order  impulse  response.  Theoreti¬ 
cally,  the  proposed  identification  technique  represents  a  signif¬ 
icant  improvement  over  existing  identification  techniques  because 
of  its  black  box  formulation.  The  intent  of  the  study  is  to  de¬ 
termine  if  this  identification  technique  can  be  practically  imple¬ 
mented  and  maintain  an  advantage  over  existing  techniques. 

The  study  effort  is  divided  into  two  parts: 

Part  I  An  implementation  feasibility  study  to  determine 

practical  methods  of  implementing  the  measurement 
scheme  -  both  digital  and  analog  -  and  to  evaluate 
the  requirements  for  the  components  of  the  measure¬ 
ment  scheme. 

Part  II  A  computational  complexity  study  of  the  identifica¬ 

tion  technique  processing  to  determine  the  class  of 
nonlinear  systems  to  which  the  technique  can  be 
practically  applied. 

This  final  report  represents  the  results  of  Part  II  of  the 
study  effort  -  the  computational  complexity  study. 

B.  SUMMARY  OF  RESULTS  AND  CONCLUSIONS 

This  part  of  the  study  effort  focused  on  identifying  the  com¬ 
putational  limitations  of  the  identification  technique  that 
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restrict  its  application  to  practical  systems  and  on  developing 
methods  of  easing  these  limitations. 

The  primary  computational  limitations  of  the  identification 
technique  arise  from  the  required  matrix  inversions  necessary  to 
evaluate  the  system  residues.  The  dynamic  range  of  the  matrix 
entries  increases  as  the  matrix  size  increases  and  these  entries 
can  violate  the  dynamic  range  constraints  of  typical  general- 
purpose  computers  even  for  moderate  size  systems.  This  problem 
is  complicated  further  when  the  linear  system  transfer  function 
is  wide  band. 


Three  approaches  are  suggested  for  alleviating  these  computa¬ 
tional  problems: 

(1)  Matrix  scaling 

(2)  Reduction  of  the  order  of  YgCs) 

(3)  Computational  scaling  for  the  total  identification  tech¬ 
nique  . 

The  third  approach  was  not  addressed  in  detail  in  the  study 
but  requires  a  continuous  scaling  of  matrix  inversion  operations 
to  take  advantage  of  the  total  dynamic  range  of  the  digital 
computer . 

Matrix  scaling  techniques  were  investigated  on  the  basis  of 
known  algorithms.  The  results  of  this  investigation  demonstrate 
that  these  techniques  offer  some  relief  of  the  computational 
problems  but  do  not  substantially  increase  the  applicability  of 
the  identification  technique. 

A  primary  conclusion  of  this  study  is  that  the  order  of  the 
second  order  response  must  be  reduced  as  much  as  possible  while 
still  maintaining  the  integrity  of  the  identification  technique. 
This  approach  dictates  a  more  extensive  measurement  process  for 
the  identification  technique.  Several  approaches  to  accomplish 
this  reduction  of  the  order  of  Y2(s)  are  postulated  in  this 
report.  These  approaches  include: 

(1)  Restricted  Frequency  Approaches  -  including  use  of  a 
low-pass  filter  at  the  system  output,  and  appropriate 
selection  of  integration  time. 


(2) 


Single  Exponential  Input  -  the  input  signal  consists  of 
a  single  exponential  function  x(t)  =  e"ait-  instead  of 


N  -a.  t 
E  e  1 
i=l 
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The  input  is  applied  N  times  (changing  a^)  and  an  appro¬ 
priate  set  of  measurements  is  taken.  The  identification 
process  is  essentially  repeated  N  times  to  generate  the 
required  set  of  linearly  independent  equations. 

(3)  Dominant  Pole  Concept  -  The  linear  transfer  function  is 
modeled  by  a  lower  order  transfer  function  where  the 
dominant  poles  are  used  in  the  transfer  function  model. 

The  computational  problem  introduced  by  a  wide-band  system, 
i.e.,  a  near  singular  matrix  to  be  inverted,  can  be  avoided  by 
the  following  method.  The  poles  of  Y2(s)  of  the  form  X*  +  Xj  =  Xj 
are  modeled  as  single  poles  at  Xj  and  the  residues  are  combined 
to  obtain  the  total  residue.  A  modification  of  the  identifica¬ 
tion  technique  allows  identification  of  the  A.,  ,  quantities. 

Klk2 

This  has  been  demonstrated  for  two  pairs  of  poles  of  the  form 


The  issue  of  the  isolation  of  the  second-order  response  from 
the  total  system  response  was  also  addressed  in  this  study.  It 
was  shown  that  it  is  not  necessary  to  isolate  y2(t)  from 
yx(t)  +  y2(t)  to  identify  the  linear  and  second-order  impulse 
responses.  However,  it  was  also  shown  that  y2(t)  must  be  isolated 
from  third  and  higher-order  system  responses  In  order  to  identify 
the  second-order  impulse  response,  hgCt^.tg). 
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SECTION  I 


INTRODUCTION 


A.  STUDY  OBJECTIVES 

The  basic  objective  of  this  study  effort  is  to  evaluate  the 
practical  feasibility  of  a  nonlinear  system  identification  tech¬ 
nique.  The  identification  procedure  studied  is  a  black  box 
technique  where  only  input  and  output  terminal  measurements  of 
the  nonlinear  system  are  used.  The  identification  technique  is 
applicable  to  a  broad  class  of  weakly  nonlinear  systems  whose 
response  can  be  characterized  by  a  finite  Volterra  series.  The 
identification  procedure  involves  processing  the  input  and  out¬ 
put  responses  of  a  nonlinear  system  to  obtain  a  set  of  linearly 
independent  equations  that  uniquely  define  the  parameters  of  a 
functional  form  of  the  second-order  impulse  response.  Theoreti¬ 
cally,  the  proposed  identification  technique  represents  a  signif¬ 
icant  improvement  over  existing  identification  techniques  because 
of  its  black  box  formulation.  The  intent  of  the  study  is  to  deter¬ 
mine  if  this  identification  technique  can  be  practically  imple¬ 
mented  and  maintain  an  advantage  over  existing  techniques. 

The  study  effort  is  divided  into  two  parts: 

Part  I  An  implementation  feasibility  study  to  determine 
practical  methods  of  implementing  the  measurement 
scheme  -  both  digital  and  analog  -  and  to  evaluate 
the  requirements  for  the  components  of  the  measure¬ 
ment  scheme. 


Part  II  A  computational  complexity  study  of  the  identifica¬ 
tion  technique  processing  to  determine  the  class  of 
nonlinear  systems  to  which  the  technique  can  be 
practically  applied. 

This  final  report  represents  the  results  of  Part  II  of  the 
study  effort  -  the  computational  complexity  study.  The  imple¬ 
mentation  feasibility  study  results  were  presented  in  Part  I  of 
this  final  report  (Reference  1). 


B.  SUMMARY  OF  RESULTS  AND  CONCLUSIONS 


This  part  of  the  study 
tational  limitations  of 


effort  focused  on  identifying  the 
the  identification  technique  that 
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restrict  its  application  to  practical  systems  and  on  developing 
methods  of  easing  these  limitations. 

The  primary  computational  limitations  of  the  identification 
technique  arise  from  the  required  matrix  inversions  necessary  to 
evaluate  the  system  residues.  The  dynamic  range  of  the  matrix 
entries  increases  as  the  matrix  size  increases  and  these  entries 
can  violate  the  dynamic  range  constraints  of  typical  general- 
purpose  computers  even  for  moderate  size  systems.  This  problem 
is  complicated  further  when  the  linear  system  transfer  function 
is  wide  band. 

Three  approaches  are  suggested  for  alleviating  these  computa¬ 
tional  problems: 

(1)  Matrix  Scaling 

(2)  Reduction  of  the  order  of  Yg(s) 

(3)  Computational  Scaling  for  the  total  identification  tech¬ 
nique. 

The  third  approach  was  not  addressed  in  detail  in  the  study 
but  requires  a  continuous  scaling  cf  matrix  inversion  operations 
to  take  advantage  of  the  total  dynamic  range  of  the  digital 
computer. 

Matrix  scaling  techniques  were  investigated  on  the  basis  of 
the  algorithms  developed  in  Reference  2.  The  results  of  this 
investigation  demonstrate  that  these  techniques  offer  some  relief 
of  the  computational  problems  but  do  not  substantially  increase 
the  applicability  of  the  identification  technique. 

A  primary  conclusion  of  this  study  is  that  the  order  of  the 
second  order  response  must  be  reduced  as  much  as  possible  while 
still  maintaining  the  integrity  of  the  identification  technique. 
In  many  instances,  this  approach  dictates  a  more  extensive 
measurement  process  than  originally  required  (Reference  1). 
Several  approaches  to  accomplish  this  reduction  of  the  order  of 
Y2  (s)  are  postulated  in  this  report.  These  approaches  include: 

(1)  Restricted  Frequency  Approaches  -  including  use  of  a 
low-pass  filter  at  the  system  output,  and  appropriate 
selection  of  integration  time 

(2)  Single  Exponential  Input  -  the  input  signal  consists  of 
a  single  exponential  function  x(t)  =  e_c‘it  instead  of 

N  -a,  t 
l  e  1 
i  =  l 
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The  input  is  applied  N  times  (changing  and  an  appro¬ 
priate  set  of  measurements  is  taken.  The  identification 
process  is  essentially  repeated  N  times  to  generate  the 
required  set  of  linearly  independent  equations. 

(3)  Dominant  Pole  Concept  -  The  linear  transfer  function  is 
modeled  by  a  lower  order  transfer  function  where  the 
dominant  poles  are  used  in  the  transfer  function  model. 

The  computational  problem  introduced  by  a  wide-band  system, 
i.e. ,  a  near  singular  matrix  to  be  inverted,  can  be  avoided  by 
the  following  method.  The  poles  of  Y2(s)  of  the  form  +  Xj  =  X 
are  modeled  as  single  poles  at  Xj  and  the  residues  are  combined  J 
to  obtain  the  total  residue.  A  modification  of  the  identifica¬ 
tion  technique  allows  identification  of  the  Aklk2  quantities. 

This  has  been  demonstrated  for  two  pairs  of  poles  of  the  form 
Xf  +  Xj  =  Xj. 

The  issue  of  the  isolation  of  the  second-order  response  from 
the  total  system  response  was  also  addressed  in  this  study.  It 
was  shown  that  it  is  not  necessary  to  isolate  y2(t)  from 
yi ( t )  +  y2(t )  to  identify  the  linear  and  second-order  impulse 
responses.  However,  it  was  also  shown  that  y2(t)  must  be  isolated 
from  third  and  higher-order  system  responses  in  order  to  identify 
the  second-order  impulse  response,  h2(t1(t2). 
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SECTION  II 


IDENTIFICATION  TECHNIQUE 


A.  IDENTIFICATION  TECHNIQUE  BACKGROUND 

The  basic  objective  of  this  study  (Part  II)  is  to  investigate 
the  computational  complexity  aspects  of  an  identification  tech¬ 
nique  for  nonlinear  systems.  The  identification  technique  is 
briefly  reviewed  in  this  section.  The  identification  technique 
is  described  in  detail  in  the  Part  I  Final  Report  (Reference  1) 
and  is  based  on  the  analysis  presented  in  Reference  3.  This 
technique  is  a  "black  box"  procedure  in  that  only  measurements  at 
the  system  input  and  output  terminals  are  required.  The  identi¬ 
fication  technique  is  applicable  to  a  class  of  weakly  nonlinear 
systems  whose  behavior  is  adequately  characterized  in  terms  of  a 
finite  Volterra  functional  series  given  by 


N  N  n 

y(t)  -  E  yn(t)  -  l  i  h  (t  . T)  n  x(t  -  T  )  dr  (1) 

n=l  n  n=l  n  1  n  p=l  P  P 


where 

yn(t)  is  the  n  order  portion  of  the  response 
1  denotes  an  n-fold  integration  from  to  +« 


n 

n  denotes  an  n-fold  product. 


p  rf  denotes  the  number  of  terms  in  the  infinite  volterra  series. 

The  n  -order  Volterra  kernel  hn(Tj_ , .  .  .  ,tn)  can  be  referred 
to  as  the  n^-order  nonlinear  impulse  response  (Reference  4),  In 


actuality,  the  nonlinear  impulse  responses  may  not  be  identically 
zero  above  order  N.  However,  the  finite  sum  of  equation  (1)  im¬ 


plies  that  higher-order  terms  contribute  negligibly  to  the  output. 


The  identification  technique  developed  in  Reference  3  is  de¬ 
signed  to  identify  the  parameters  of  closed-form  expressions  for 
the  nonlinear  impulse  responses,  h  (tj., 
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n  ■  1,  2,...,N.  The  analysis  presented  in  Reference  3  demon¬ 
strates  how  the  technique  identifies  the  parameters  of  hi(t), 
h2(ti,  t£)  and  t^tj,  t2»  t2).  On  the  basis  of  this  analysis,  it 
is  believed  that  the  technique  is  extendable  to  identification  of 
higher  order  nonlinear  impulse  responses  (N  >_  4).  This  study 
(Part  II)  is  concerned  with  the  computational  aspects  of  the  iden¬ 
tification  of  the  linear  and  second-order  nonlinear  impulse 
responses,  hi(t),  h2(ti,t2)- 

It  has  been  shown  (Reference  3)  that,  when  the  linear  incre¬ 
mental  model  of  a  nonlinear  system  is  described  by 


hx(t )  « 


N  At 

E  R.e  1 
i=l  1 


,  t  >  0 


,  t  <  0 


(2) 


where  Re  { X ^ }  <  0  and  it  is  assumed  that  the  A^  are  distinct, 
the  second-order  nonlinear  impulse  response  can  be  expressed  in 
the  symmetrical  form  (Reference  3): 


M  N 

h2(ti,t2)  *  £  I  Ak,k0  e 

kj=l  k2=l  1  2 


Vtl+akQt2 

1  z 


v(tn  - 


M  N  ak1t2+ak9tl 

+  l  E  Ak  .  e  x  *  U(tn  -  t9) 
ki-l  k2-i  kik2  1  2 


(3) 


where 


M  ®  N  +  1. 


(4) 


U(t)  = 


,  t  >  0 
,  t  <  0 


(5) 


and  where  the  natural  frequencies  in  equation  (3)  are  related  to 
those  in  equation  (2)  according  to: 


al  *  Xl»  a2  =  X2 . aN  *  XN 

aN+l  =  X1  "  A1  =  °’  aN+2  =  \±  -  A2,...,  a2N 

a2N+l  X2  Xl’  a2N+2  "  X2  “  X3*-'**  a3N-l 
aN2-N+3  XN  Xl*  aN2-N+4  =  XN  X2,,••, 
aN2+l  XN  XN-1‘ 


The  ordering  of  the  terms  in  equation  (3)  assumes  all 
the  factors  A^  -  A jto  be  distinct,  such  that  Ai  -  A j  f  A^  for  any 

i,j,k  =  1 . N.  Also,  the  zero  entry  that  results  from  A^  -  Aj 

when  i  =  J  is  included  only  once  as  the  entry  aN+i-  In  addition, 
it  is  readily  shown  that  (Reference  3) 


klk2  ^2^1 


for  <  N 


(V) 


and  that  the  coefficients  of  terms  in  equation(3)  having  the  form 


(A 

e 


i 


V‘l  +  V2 


i  f  1 


are  identically  zero. 

The  identification  technique  identifies  the  parameters  of 
h2(ti,  t2>  as  represented  in  equation  (3). 

B.  IDENTIFICATION  TECHNIQUE  DESCRIPTION 

The  functional  form  for  t^Ct^,  t2)  established  in  equation 
(3)  implies  that  the  identification  of  h2(tj,  t2>  reduces  to  iden¬ 
tification  of  the  parameters  a^  ,  aj~  ,  A^,^  and  N.  However, 
equations  (4)  and  (6)  show  that  ak.  ,  ak  ana  N  can  be  determined 
once  the  linear  impulse  response  is  known.  Therefore,  the  task 
of  identifying  these  parameters  reduces  to  the  task  of  identify¬ 
ing  hjCt).  The  problem  of  identifying  the  coefficients  A^ 
still  remains.  1  ^ 

The  identification  process  separates  into  two  distinct  steps: 
(1)  identification  of  h^ ( t ) ;  and  (2)  identification  of  the  A^j^ 
quantities  of  h2<ti,  t2)-  These  two  steps  are  considered  below. 

1.  Identification  of  the  Linear  Impulse  Response,  h]_(t). 

The  first  step  in  the  identification  of  hj.(t),  the  linear 
impulse  response  of  a  nonlinear  system,  is  to  excite  the  system 


with  an  input  amplitude  such  that  the  output  is  linear.  The  amp¬ 
litude  of  this  signal  can  be  determined  by  exciting  the  system 
with  a  sinusoidal  signal  of  amplitude  A  and  performing  a  spectral 
analysis  of  the  resultant  response.  Amplitude  A  is  then  adjusted 
until  the  amplitude  level  of  the  harmonic  frequencies  of  the 
output  becomes  sufficiently  small  compared  to  the  level  of  the 
fundamental  component.  The  poles  and  residues  of  hi(t)  be  will 
identified  using  the  pencil-of-functions  approach  (Reference  5). 
The  pencil-of-functions  approach  integrates  the  input  to  the 
linear  system  and  resulting  output  N  times  over  the  real-time 
interval  (0,T). 


It  has  been  shown  (Reference  5)  that  poles  of  the  linear 
system  satisfy  the  polynomial  equation 


(8) 


where  G2JJ+1  is  the  Gram  determinant  shown  in  equation  (9)  below: 


2N+1 


A 

H 

I  >> 

H 

l>> 

V 

<yl'y2> 

<yi’iw 

<yl'x2>  • 

•  <yl’XN+l 

'W  ••• 

<y2,yN+l> 

<y2'x2>  • 

•  <y2  ’  XN+1 

<yN'V 

<yN’V 

<yN ' yN+l> 

<yN’V  • 

<yN ’ XN+1 

<x2 ’ yl> 

<x2’y2>  ••• 

<x2,yNH  > 

<x2’ V  1 

••  <X2 ’ XN+1 

<xN+l’yl><XN+l’y2>- ’ • <xN+l,yN+l><xN+l,x2>' ' ,<XN+1,XN+1> 

(9) 


and  where 


xl+l(t)  = 


/x. (x )dt 
0 


0  <  t  <_  T 
elsewhere 


(10) 


i  =  1 . N 


7 


I  y,(T)dT 


y^i  (t) 


0  <  t  <  T 


elsewhere 


(11) 


Further,  the  residues  R^  of  the  poles  Xi  satisfy  the  equation 


R  «  C'1Y 


(12) 


where 


R  =  residue  matrix  =  Ir, 


(13) 


y2(T) 

»3<T> 

Y  =  output  matrix  ■  y4(T) 


(14) 


C  =  N  x  N  matrix  whose  i,j  element  is  defined  by 


P1<T>  i  WT> 

Tf  '  m=l  (X  >1+i- 


(15) 


where 


T  A,(T-t) 

P  .  (T)  =  /  e  J  x(x )dx 

•J  n 


(IB) 
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References  1  and  3  describe  how  this  processing  can  be  used  to 
to  determine  N. 


2.  Identification  of  the  Second  Order  Impulse 
Response, 

The  second  step  of  the  identification  procedure  is  to 
identify  the  unknown  parameters  of  hgC^,  tg).  hj^t^.tg) 


given  by: 


M 


N  VV ftk„t2 

h„(t1  ,t  )  =  Z  Z  A.  e  A 

*  1  *  k1=l  k2=l  k1k2  ~~ 


ak  Vak  ti 

k  k  e  1  2  u(ti  -  V 

kj-1  k2-l  kik2  1  * 


M 

+  1 


N 

Z  A, 


(17) 


the  only  unknown  parameters  are  the  A^j^  quantities  since  M,  N, 
Afci  and  A^q  are  known  from  identification  of  hi(t).  A  procedure 
for  determining  the  Afc.j^  using  the  pencil-of-functions  method  is 
described  in  this  section. 


The  identification  procedure  utilizes  the  response  of 
the  weakly  nonlinear  system  to  a  sum  of  L  decaying  exponentials 
as  described  by: 


x(t) 


L  -a .  t 
Z  e  1 
i=l 


t  >  0 

t  <  0 


(18) 


where  Re  { a ^ )  >  0.  The  second-order  portion  of  the  response  to 
x(t)  is  given  by 


M  N  L  L 


where 


The  expression  in  equation  (19)  is  the  Laplace  trans¬ 
form  of  a  cum  of  exponential  time  functions.  This  sum  can  be 
interpreted  as  the  impulse  response  of  an  equivalent  linear  sys¬ 
tem  as  indicated  in  Figure  1.  In  other  words,  the  second-order 
response  y2(t)  can  be  visualized  as  though  it  were  generated  by 
an  equivalent  linear  system.  However,  the  equivalence  is  valid 
only  if  the  equivalent  linear  system  is  considered  to  be  excited 
by  an  impulse.  It  follows  that  the  problem  of  identifying 
h2(tl>t2^  has  been  reduced  to  the  simpler  problem  of  identifying 
a  linear  system  and  the  pencil-of-functions  technique  can  be  used 
again . 
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I 


.1 


{y2(»)} 


Figure  1.  Equivalent  Linear  System  with  Transfer 

Function  Yg(s) 


The  system  is  excited  by  an  input  amplitude  such  that 
the  output  is  described  by  linear  and  second-order  terms,  y^(t) 
and  y2(t ).  The  identification  process  will  operate  on  the  signal 
y2(t).  For  this  purpose,  the  second-order  portion  of  the  re¬ 
sponse,  y2(t),  is  isolated  from  the  total  response.  y2(t)  is 
obtained  by  subtracting  from  the  total  response,  the  correspond¬ 
ing  linear  response  yi(t),  which  is  known  because  hi(t)  has  been 
identified.  It  is  shown  in  Reference  2  and  in  Section  III.F  of 
this  report  that  the  second-order  response  y2<t)  need  not  be  iso¬ 
lated  from  the  total  response  for  the  identification  procedure  to 
work.  However,  isolation  of  y£(t)  from  the  total  response  eases 
the  mathematical  presentation  and  is  assumed  at  this  point. 

Once  y2(t)  is  isolated  from  the  total  response,  the  co¬ 
efficients  are  then  evaluated  by  applying  the  pencil-of- 

functions  method  to  y2(t),  treating  it  as  though  it  were  the  im¬ 
pulse  response  of  a  linear  system.  This  latter  step  is  now  dis¬ 
cussed  in  detail. 

From  equation  (19),  the  poles  of  Yg(s)  are  given  by 
s  =  a,  +  a.  ,  k-  =  1,  .  .  .  ,M;  kQ  =  1 . N 

kl  k2  1  * 

s  =  -a.  +  a.  ,  i  =  1,...,L;  k0  =  1,...,N 


s  =  -ai  -  otj  ,  i,j  *  1,...,L. 


(21) 
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•*  I 

I 


First,  consider  poles  of  the  form  s  *  a^-  +  a^  =  2  A*;  £ 

•  The  terms  in  Y2(s)  corresponding  to  the  pole  at  2  A^ 

are  given  by 


L  L 


ws>  - *  * A 


w  ji i A**  <«j *  v<°i  * 


If  the  residue  of  the  pole  at  2  X£,  as  evaluated  using  the  pencil- 
of-functions  method,  is  it  follows  that 


L  L  1 


l  *  1,...,N. 


This  procedure  results  in  identification  of  N  of  the  coefficients. 

Consider  next  poles  of  the  form  s  =  a^i  +  aic2  =  A £  +  A 
where  £  f  m  and  £,m  =  Since  A£m  =  A^  for  £ ,  m  <_  N, 

the  terms  in  Y2(s)  corresponding  to  the  pole  at  A^  +  Aro  are  given 
by 

Y2£m(s)  =  =  =  A£m 

i=l  j=l 

+  2A^ 

(ccj  +  X£)(ai  +X(!)(ai  +  +  A£  +  Am) 

Of  +  a .  +  2A  ^  1 

+  <“j  +  +  V<at  ”aJ  *  X1  '  h  -  Xn)J  •  <24> 

If  the  residue  of  the  pole  at  A^  + A  m,  as  evaluated  using  the 
pencil-of-functions  method,  is  e ^  ,  it  follows  that 


A£m  ~  ^£m 


l  l  r 

li-1  0=1  (aj  + 


ai  +  ai  +  2X£ 

A  £ )  (a  j  +  V(ai  +  aj  +  X£  +  V 


_ ai  +  aj  +  2Xm _ 

(aj  +  Xm)(ai  +  Xm)(ai  +  aj  +  X£  +  Xm>  £,m  =  1 . N 

•if  m ,  £  <  m. 


namrms^ 


i*j*  *  t  ■ 


This  procedure  results  in  identification  of  N(N  -  l)/2  of  the 
coefficients. 

o 

The  remaining  unknown  N  coefficients  cannot  be  evalu¬ 
ated  directly,  as  was  done  in  equations  (23)  and (25),  because 
the  residues  of  the  other  poles  in  Yo(s)  involve  linear  combina¬ 
tions  of  more  than  one  unknown  coefficient.  However,  if  the 
number  of  exponential  input  signals,  L,  is  set  equal  to  N,  N2 
linearly  independent  equations  involving  the  N2  unknown  A^llC2 
coefficients  can  be  obtained  by  considering  the  poles  of 
Y2(s)  of  the  form  s  »-Oj  +  i  85  1,...,N;  j  *  1,...,N,  in  a 
manner  similar  to  the  above  analysis.  This  fact  is  proven  in 
Reference  3.  Solution  of  the  N2  equations  completes  the  identi¬ 
fication  process. 
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SECTION  III 


SYSTEM  COMPLEXITY  STUDY 


A.  COMPUTATIONAL  COMPLEXITY  CONSIDERATIONS 

The  analyses  of  Section  II  and  Reference  3  demonstrated  that 
the  theoretical  derivation  of  the  identification  technique  was 
restricted  to  a  class  of  nonlinear  systems  described  by  the  im¬ 
pulse  response  of  the  form 

N  X.t 

Z  R.  e  1  ,  t  >  0 
i=l 

h1(t)  =  (26) 

0  ,  t  <  0 

Although  there  appears  to  be  no  theoretical  limitation  preventing 
the  application  of  the  identification  technique  to  systems  with 
multiple-order  poles,  the  analysis  has  not  been  done  to  support 
this  conclusion.  The  systems  modeled  by  equation  (26)  represent 
a  broad  class  of  nonlinear  systems  to  which  the  technique  can  be 
applied.  Practical  limitations  of  the  technique  will  restrict 
the  class  of  systems  to  a  subset  of  those  represented  by  equa¬ 
tion  (26).  These  practical  limitations  arise  primarily  from  the 
computational  requirements  of  the  identification  technique  proc¬ 
essing  scheme.  These  limitations  constrain  the  maximum  value  of 
N,  which  restricts  application  of  the  technique  to  systems  whose 
linear  incremental  model  has  N  poles  or  less. 

This  section  investigates  these  computational  limitations, 
attempts  to  establish  a  maximum  value  for  N,  and  presents 
selected  techniques  to  alleviate  these  computational  problems. 

1.  Computational  Complexity  Limitations 

The  numerical  computation  requirements  of  the  identi¬ 
fication  technique  are  summarized  below.  They  are: 

(1)  N  numerical  integrations  of  input  and  output 
(N  is  the  order  of  the  system). 

(2)  Formation  of  2N  +  1  inner  product  entries  for  the 
Gram  matrix. 
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(3)  Evaluation  of  determinants  of  N(2N  +  1)  matrices. 

t*  h 

(4)  Solution  of  a  order  polynomial  equation. 

2 

(5)  Evaluation  of  the  N  entries  of  the  C  matrix 

in  the  residue  equation. 

(6)  Inversion  of  an  N  dimension  C  matrix. 

(7)  Solution  of  the  quantities  (second-order 

system) . 

The  numerical  accuracy  requirements  for  these  computa¬ 
tions  were  investigated  in  Part  I  of  this  study  (Reference  1). 

The  numerical  accuracy  required  for  satisfactory  performance  of 
the  identification  technique  increases  significantly  with  in¬ 
creasing  N.  The  severest  computation  requirements  are  imposed  by 
(1)  the  formation  and  inversion  of  the  C  matrix  used  in  the 
residue  equation  R  =  C”1  Y;  and  (2)  numerical  integration  and 
formation  of  the  inner  products  for  the  appropriate  Gram  matrix. 
These  two  areas  are  addressed  below. 

2.  Residue  Equation  Computational  Requirements 

The  residue  equation  for  the  identification  technique 
is  given  by 

R  =  C'1  Y  (27) 


The  C  matrix  has  dimension  N'  where  N'  is  the  order  of  the  sys¬ 
tem  being  identified.  The  C  matrix  entries  are  given  by 


j 


T  -X,x 

f  e  J  x(x)  dx 
0 


(28) 


where  the  X,,  j  =  1 . N'  are  the  poles  of  the  system  being 

identified. 


or 


For  the  identification  technique,  x(t)  is  of  the  form 
-a.t 

x(t)  =  e  ,  t  >  0  (linear  system  identif ication)(29) 
x(t)  =  6 ( t )  (secondary  system  identification)  (30) 
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(2)  for  x(t)  *  6(t) 


(32) 

It  has  been  noted  in  previous  work  (References  2,  4) 
that  the  C  matrix  tends  to  be  ill-conditioned,  which  hampers  its 
computational  inversion.  Two  significant  problems  complicate 
the  inversion  problem:  (1)  the  dynamic  range  of  the  C^j ;  and 
(2)  the  near  singularity  of  the  C  matrix  when  two  of  the  N'  poles 
are  nearly  equal.  The  singularity  problem  is  addressed  in  Sec¬ 
tions  III.B  and  C  while  the  dynamic  range  problem  is  addressed 
below. 

The  computational  requirements  of  inverting  the  C 
matrix  are  basically  determined  by  its  dimension.  The  dimen¬ 
sion  of  the  C  matrix  is  determined  by  the  number  of  poles  of  the 
system  being  identified.  Consider  a  nonlinear  system  whose 
linear  transfer  function  has  N  poles.  The  dimension  of  the  C 
matrix  is  then N  x  N.  Identification  of  the  second-order  trans¬ 
fer  function  involves  the  identification  of  the  residues  of  the 
second-order  response,  Y2(s).  The  number  of  poles  of  Y2(s)  is 

N'  *  N  (|  +  |  +  L)  +  L  ^  -2  ^  (33) 

where  L  is  the  number  of  exponential  signals  composing  the  input. 
For  the  general  case  where  L  =  N, 

N*  ■  2N  (N  +  1)  (34) 
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Figure  2  plots  N'  as  a  function  of  N.  It  is  noted  that,  for  a 
10-pole  linear  system,  the  calculation  of  the  residues  of  Y2(s) 
involves  inversion  of  a  240  x  240  matrix.  The  dimension  of  the 
C  matrix  for  second-order  identification  grows  rapidly  with  the 
number  of  linear  system  poles. 


The  dynamic  range  of  the  C^.,  entries  becomes  signifi¬ 
cant  as  the  dimension  of  the  C  matrix1-1 increases  since  the  Cij 
entries  are  inversely  proportional  to  Xj1 ,  for  i  ■  This 

is  demonstrated  below. 


For  identification  of  the  residues  of  Y2(s),  the  input 
x(t)  is  given  by 

x(t)  =  Ra  6 ( t )  (35; 

where  6(t)  is  the  unit  impulse.  The  Cij  entries  for  this  input 
are  given  by  r  .  _  n 

|e  J  i  x®  -  1 

J  [Xj1  m=l  (Xj)1  1  '  m  (m  -  1) ! J 

Suppose  the  system  response  of  interest.  Y2(s),  has  a  maximum 
pole/minimum  pole  ratio  -  10.  Further,  assume  that 


or  that 


-10  <  XjT  <  -1 

T  -  -1/Xmin. 


For  the  minimum  pole  Xj  min 


Ra  r — r  -  r  — ^ — l 

a  L  (-1  y  m-l  (-1)1  (in  -  1)  !  J 


(37) 


For  different  values  of  i,  is  given  by 


1.37  R 


2  -1.63  R, 


10  -0.35  R, 


N.  NUMBER  OF  UNEAR  SYSTEM  POLES 


Figure  2.  Number  of  Poles  of  Y2(s)  as  a  Function  of  the 
Number  of  Linear  System  Poles 
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For  the  maximum  pole  x 


ij 


j  max 


Te'10  1 

- - 7  -  E  — 

.(-10)1  m=l  (- 


IzlI 


m  -  1 


10)1  +  1  -  m  (m 


-J 


(38) 


Again,  for  different  values  of  i,  is  given  by 


i 

Cij 

1 

-0.1 

2 

0.09 

10 

io-8 

io 

io-10 

The  dynamic  range  of  the  entries  for  a  C  matrix  of 

dimension  100  is  approximately  lO^OO .  For  N'  greater  than  100, 
the  dynamic  range  is  even  greater.  This  dynamic  range  can  cause 
significant  difficulty  when  the  inverse  of  the  C  matrix  is  evalu¬ 
ated.  For  a  typical  general-purpose  computer,  the  maximum  com¬ 
putable  dynamic  range  is  IO7*5  (10“3  to  10  ).  Furthermore, 

matrix  inversion  involves  multiplication  and  division  of  pairs  of 
matrix  entries.  The  resultant  product  or  division  must  be  in  the 
allowable  dynamic  range,  which  implies  that  the  individual  matrix 
entries  must  be  well  within  the  dynamic  range. 

For  a  system  with  two  poles  of  its  linear  transfer 
function  with  a  ratio  of  10,  the  computer  limitation  constrains 
the  class  of  systems  to  those  with  N  <  4.  This  problem  can  be 
alleviated  somewhat  with  scaling  but  not  significantly. 

A  method  of  reducing  the  dimension  of  the  C  matrix  must 
be  found  to  ease  the  computational  limitations.  This  requires 
that  the  number  of  poles  Y2(s)  be  reduced. 

A  significant  reduction  in  the  number  of  poles  of  Y2(s) 
results  if  L  =  1  instead  of  L  =  N,  where  L  is  the  number  of  ex¬ 
ponentials  used  in  the  input  function.  For  L  =  1,  the  number  of 
poles  of  Y2(s)  is  given  by 

f  i  (39) 


N'  =  N 


,  N  .  5. 
(  2  + 
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This  is  plotted  in  Figure  3  along  with  the  plot  for  L  =  N.  It  is 
noted  that  the  number  of  poles  of  Y2(s)  is  significantly  less  for 
L  -  l . 


It  is  recalled  that  L  »  N  was  required  in  order  to 
generate  a  linearly  independent  set  of  equations  from  the  resi¬ 
due  of  the  poles  at  s  =  +  a^o ,  i  =  1, . . . ,L:  k£  11  1,  .  .  .  ,N. 

These  equations  can  also  be  obtained  by  using  L  «  1  and  exciting 
the  system  with  N  individual  inputs  and  recording  the  response 
to  each  input.  This  obviously  complicates  the  identification 
procedure  but  does  reduce  the  magnitude  of  the  matrix  inversion 
problem. 

If  the  approach  is  adopted  (L  =  1),  then  the  identifi¬ 
cation  procedure  is  modified  as  follows.  The  system  is  excited 

by  the  input  x(t)  =  e  alt.  The  resultant  output  is  of  the  form 


M  N  r _ 2 

1»1  k2=l  klk2  [(<xl  +  +  ak1 

(s  -  V*  v) 


Y9(s)  =  E 
d  k 


+  a.  ) 
k2 


'  («!  ♦  »klX«i  *  + 

+  _  2 _  1  1 

<°i  *  ak2)(2“i  +  akj  +  ak2>(s  *  2aiJ 


(40) 


The  number  of  poles  in  Y2(s)  above  is  (N/2)(N  +  5)  +  1.  The 
Aklko  quantities  of  the  form  A^,  i  =  are  identified 

directly  from  the  residues  of  tne  poles  at  s  =  2ai.  The 
Akik2  quantities  of  the  form,  A^m,  S-  t  m,  l,  m  *  1,...,N  are 
identified  directly  from  the  residues  of  the  poles  at  s i  =  a£  +  am- 
This  procedure  identifies  N  +  [(N  -  1)/2]N  Aj^i^  quantities. 

-otgf 

The  system  is  then  excited  by  the  input  x(t)  =  e 
The  resultant  Y2(s)  has  (N/2)(N  +  5)  +  1  poles  but  the  contribu¬ 
tions  from  N  +  £N  -  l)N/2]  poles  are  known  and  may  be  subtracted 
out  from  the  second-order  response.  The  resultant  Y2(s)  has 
2N  +  1  poles.  The  resultant  C  matrix  is  reduced  to  order  2N  +  1. 
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Figure  3.  Number  of  Poles  of  Y2(s)  as  a  Function  of  t 
Number  of  Linear  System  Poles  and  Exponential  Inputs 


-a^t 

This  process  is  repeated  for  each  input  x(t)  *  e  , 
i  =  3,...,N.  Each  process  results  in  a  set  of  equations  involv¬ 
ing  the  residue  at  the  poles  s  =  -  +  akg*  The  first  applica¬ 

tion  involves  inverting  an  (N/2)(N  +  5)  +  I  matrix  while  the  re¬ 
maining  (N  -  1)  applications  require  inversion  of  a  2N  +  1  matrix 

The  advantage  of  this  approach  is  the  reduction 
achieved  in  the  number  of  poles  of  Y2(s).  This  eases  the  compu¬ 
tational  complexity  of  matrix  inversion.  The  primary  disadvant¬ 
age  of  the  technique  is  that  the  identification  measurement  proc¬ 
ess  must  be  repeated  N  times  in  order  to  identify  l^Ct^tg). 

Another  alternative  approach  to  reducing  the  number  of 
poles  of  Y2(s)  without  significantly  increasing  the  measurement 
process  is  to  initially  excite  the  system  with  an  input,  x(t)  31 

e-0!*.  The  residues  of  the  poles  at  s  *  +  Xj ,  i,  j  = 

are  evaluated  to  determine  the  appropriate  .  quantities.  The 
system  is  then  excited  with  an  input  K1K2 

N  -a .  t 
x(t)  =  E  e 
i  =  l 


Since  N  +  [N(N  -  l)/2]  Ak^k2  quantities  have  been  identified 
above,  the  contributions  of  the  associated  poles  can  be  sub¬ 
tracted  from  the  total  response.  The  resultant  response  has 
( 3N/2 )  (N  +  1)  poles.  This  is  a  reduction  of  25  percent  or 
(N/2)  (N  +  1)  from  the  original  second-order  response  obtained 
with 


N  -a ,  t 
x(t)  =  I  e 
i=l 

Although  the  order  of  reduction  is  not  as  great  as  that 

achieved  by  applying  the  input  x(t)  =  e'"ait'  N  times,  the  identi¬ 
fication  measurement  process  need  be  repeated  only  twice. 

The  basics  of  these  approaches  to  the  identification 
process  are  summarized  in  Table  1  for  comparison  purposes. 

The  best  approach  for  the  identification  technique  is 
dependent  on  the  system  under  test.  The  number  of  linear  system 
poles  and  the  ratio  of  the  maximum  to  minimum  poles  dictate  the 
complexity  of  the  matrix  inversion  problem  and,  in  turn,  deter¬ 
mine  which  of  the  above  approaches  will  maximize  the  performance 
of  the  identification  technique.  Therefore,  selection  of  the 
best  approach  must  wait  until  the  linear  portion  of  the  system 
under  test  has  been  identified. 
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Another  potential  approach  to  reducing  the  computa¬ 
tional  complexity  is  the  frequency  range  separation  technique  of 
Jain  and  Osman  (Reference  2).  This  approach  divides  the  fre¬ 
quency  extent  of  a  given  system  into  a  set  of  frequency  ranges, 
e.g.,  low  frequency  region,  middle  to  high  frequency  transition 
region,  and  high  frequency  region.  The  system  is  excited  by  an 
input  signal  that  is  approximately  matched  to  the  frequency 
region  of  interest  and  the  integration  time  is  selected  consist¬ 
ent  with  this  frequency  range.  The  identified  transfer  function 
is  then  a  representation  of  the  system  transfer  function  in  the 
specified  frequency  region. 

This  approach  assumes  some  a  priori  knowledge  of  the 
poles  of  the  system  transfer  function  in  order  to  permit  fre¬ 
quency  region  separation  and  determination  of  the  number  of  poles 
of  the  system  transfer  function  in  each  region.  Since,  for  the 
nonlinear  system  identification  technique  of  interest  in  this 
study,  the  order  and  pole  locations  of  the  second-order  response 
are  known,  the  above  requirement  is  satisfied.  The  input  func¬ 
tion  for  the  nonlinear  system  identification  technique, 

N  -a .  t 
x(t)  =  E  e  1  ’ 
i=l 


could  be  divided  into  the  frequency  regions  of  interest  and  ap¬ 
plied  separately  for  each  frequency  region.  For  example,  if 
N  =  6,  and  Y2(s)  is  divided  into  three  frequency  ranges,  then 
the  identification  procedure  is  conducted  as  follows.  Three 
sets  of  measurements  are  taken,  each  with  input 

2  -a  .  t 
x(t)  =  E  e 
i=l 

where  the  ai  are  selected  consistent  with  the  frequency  region 
of  interest.  The  three  sets  of  measurements  are  then  collec¬ 
tively  used  to  solve  for  the  Akj_k2  quantities  in  the  normal 
manner. 


The  achievable  reduction  in  computational  complexity 
using  this  approach  is  dependent  on  the  characteristics  of  the 
system  under  test.  If  the  poles  of  Y2(s)  are  distributed  uni¬ 
formly  in  frequency,  then  the  number  of  poles  of  Y2(s)  in  each 
of  three  frequency  regions  is  N'/3  and  a  one-third  reduction 
in  the  size  of  the  matrices  to  be  inverted  has  been  realized. 

It  is  necessary  to  point  out  that  the  price  of  this  reduction 
is  the  need  to  perform  the  identification  procedure  three  times 
instead  of  once. 
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3. 


Numerical  Integration  Techniques 


The  results  of  the  implementation  feasibility  portion 
(Part  I)  of  this  study  indicated  that  significant  accuracy  was 
required  of  the  numerical  integration  technique  to  achieve  sat¬ 
isfactory  technique  performance.  The  results  indicated  that  the 
primary  source  of  inaccuracy  in  the  integration  was  the  quanti¬ 
zation  error  introduced  by  the  A/D  converter.  This  will  be  the 
case  independent  of  the  numerical  integration  technique  used. 
However,  the  results  of  Part  I  of  the  study  also  indicated  that 
the  Simpson's  rule  of  numerical  integration  technique  introduced 
numerical  inaccuracy  for  higher-order  systems.  The  reason  for 
this  inaccuaracy  was  the  fact  that  Simpson's  rule  reduces  the 
number  of  samples  on  each  successive  integration.  Simpson's 
rule  of  integration  (Reference  2)  is  given  by 


b 

/  y(t)  dt  = 
a 


(b  -  a) 
6  n 


y(o> 


+  ...  +  2y( (2n 


+  4y(A  T)  +  2y(2AT)  +  4y(3AT) 
-  2 ) AT )  +  4y((2n  -  1)AT) 


+  y(2nAT) 


(41) 


where 


AT  =  (b  -  a)/2n  *  time  between  samples 

2n  =  number  of  subintervals  between  data  points. 

Each  integration  using  the  Simpson's  rule  integration 
technique  results  in  a  reduction  of  the  number  of  samples  that 
can  be  used  for  the  next  integration.  This  is  illustrated  below. 

Consider  the  output  samples  yi(0),  yi(T),  yi(2T),..., 
y^(2nT),  where  nT  is  the  n*h  sample  and  T  is  the  sampling  inter¬ 
val.  The  integral  of  yi(t),y2(t),  as  obtained  using  Simpson's 
rule,  is  given  by  the  samples 

y2(°).  y2(2T),  y2(4T) , . . . ,y2(2nT) . 

It  is  noted  that  there  are  only  nT  samples  of  y2(t) 
whereas  there  were  2nT  samples  of  yj(t).  As  this  output  is  suc¬ 
cessively  integrated,  the  time  distance  between  samples  increases 
and  the  numerical  accuracy  of  the  integration  technique  decreases 
This  will  have  an  adverse  effect  on  the  performance  of  the  ident¬ 
ification  technique,  especially  for  higher-order  systems. 
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A  way  of  alleviating  this  decreasing-number-of-samples 
problem  is  to  interpolate  between  samples  output  from  Simpson's 
rule  of  integration.  For  instance,  in  the  above  example,  Simp¬ 
son's  rule  produced  an  output  at  t  *  0  and  t  *  2T,  namely, 
y2(0)  and  y2(2T).  Interpolating  linearly  between  samples  yields 


y2m 


y2(2T)  +  y2(0) 
2 


(42) 


To  determine  the  impact  of  this  procedure  on  the  performance  of 
the  identification  technique,  this  procedure  was  added  to  the 
computer  simulation  of  the  identification  technique.  The  com¬ 
puter  simulation  of  the  identification  technique  was  discussed 
in  detail  in  Part  I  of  this  final  report  (Reference  1).  The  sim¬ 
ulation  was  run  for  two  systems,  a  two-pole  and  a  four-pole  sys¬ 
tem,  The  results  are  presented  in  Tables  2  and  3.  The  perform¬ 
ance  of  the  direct  application  of  Simpson's  rule  is  included  for 
comparison  in  Tables  2  and  3, 

The  results  of  Table  2  indicate  that,  for  a  two-pole 
system,  the  interpolation  procedure  slightly  degrades  the  per¬ 
formance  of  the  technique.  The  results  of  Table  3  indicate  that, 
for  a  four-pole  system,  the  interpolation  scheme  offers  slightly 
improved  performance  for  A/D  converters  with  20  bits  or  less  of 
resolution.  This  improvement  will  continue  to  be  evident  as  sys¬ 
tem  order  increases.  These  results  suggest  that  the  numerical 
integration  technique  be  modified  to  include  this  interpolation 
method. 


B.  MATRIX  INVERSION/SCALING  TECHNIQUES 


The  primary  computational  problem  of  the  nonlinear  system 
identification  technique  is  the  matrix  inversion  involved  in 
solving  the  residue  equation. 


The  matrix  to  be  inverted  has  entries  given  in  general  by 


where 


'ij 


Pj(T) 


'm  +  1 


(T) 


1  "  m=l  + 


1  -  m 


T  X  ,(T  -t) 

P,  (T)  =  1  e  J  x ( t )  dx 

J  0 


(43) 


(44) 


For  the  second-order  system  identification,  x(t)  =  6(t),  a  unit 
impulse  which  reduces  the  matrix  entries  to 
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TABLE  3.  IMPACT  OF  INTERPOLATION  BETWEEN  SAMPLES  ON  INTEGRATION  FOR  FOUR  POLE  SYSTEM. 

INTEGRATION  TIME  =  4.8  |iS,  SAMPLING  TIME  =  1.5  ns 
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C, 


,V  i 


r,m  -  1 


lJ  X  j1  m“1  (m  -  1) ! ^  +  1  ■  m) 
Thts  expression  can  be  rewritten  as 


'ij 


1  _ 


XjT  i  ( X  jT) 


m  -  1 


But  since 


V  __  -  OjT) 


'  m»l  <»  *  D' 


m  -  1 


=  I 
m=l 


(m  -  1)! 


for  Xj  real,  the  entries  become 


'ij  i 

Xj 


«  ( X ,T) 

Z  — J 

m=i+l 


m  -  l-l 
(m  -  1)! 


or 


,  •  (XjT) 

c  =  -  -  y  -J- 

ij  x  i  tn! 

A  j  m-  i 


m 


(46) 


(47) 


(48) 


(49) 


This  expression  serves  to  illustrate  two  basic  problems  with  the 
numerical  inversion  of  the  C  matrix.  First,  the  dynamic  range 
limitation  of  a  digital  computer  limits  the  number  of  terms  which 
can  be  used  in  a  given  summation.  Also,  it  is  clear  that  if  two 
Xj  quantities  used  in  the  Cjj  expression  are  nearly  equal,  the 
Cjj  terms  become  nearly  equal  since  the  m!  quantity  tends  to  re¬ 
duce  the  difference  between  the  C^j  entries. 


1  This  numerical  similarity  between  the  Cij  entries  causes  the 

major  problem  incurred  when  attempting  to  invert  the  C  matrix. 

*  Several  approaches  to  alleviating  thisproblem  have  been  postulated 
(Ref.  2)  and  these  provide  some  relief  but  do  not  eliminate  the 

*  problem.  The  dynamic  range  of  the  matrix  entries  provides  the 
major  limitation  and  avoiding  this  problem  requires  intricate 
programming  which  essentially  corresponds  to  scaling  quantities 
after  each  operation.  This  is  a  long  and  involved  process,  and 
the  design  of  such  a  program  is  beyond  the  scope  of  the  present 
effort.  However,  it  remains  an  important  area  to  consider  in 
future  efforts  because  it  may  help  to  eliminate  a  major  limita¬ 
tion  of  the  application  of  the  technique.  The  other  techniques 
discussed  in  Reference  2  have  been  used  to  invert  matrices  which 
have  the  numerical  structure  of  the  C  matrix  and  they  offer  some 
relief  to  the  existing  problem.  These  techniques  are  reviewed 
below. 
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First,  it  should  be  noted  that  standard  matrix  inversion 
techniques  perform  well  when  the  C  matrix  is  not  nearly  singular 
and/or  large  enough  that  the  dynamic  range  limitations  of  the 
computer  are  exceeded.  This  was  demonstrated  in  References  1 
and  3  for  two-,  four-  and  eight-pole  systems.  The  postulated 
approaches  should  be  applied  only  when  the  standard  techniques 
fail  to  perform  satisfactorily. 


The  initial  technique  to  aid  in  matrix  inversion  is  called 
"adaptive  scaling"  (Reference  2).  This  technique  depends  on 
row  and  column  scaling  to  alleviate  the  problem  caused  by  matrix 
entries  which  differ  widely  in  magnitude.  If  it  is  desired  to 
invert  a  matrix  C,  the  first  step  is  to  do  a  row  and  column 
scaling  on  C,  transforming  it  to 

’’  C0  «=  PCQ  (50) 

where  P  and  Q  are  diagonal  scaling  matrices.  The  entries  of  the 
P  and  Q  matrices  can  be  chosen  as  follows.  Consider  the  P 
matrix.  The  Pj^  entry  is  computed  as  follows: 

Plt  -  n  {<c)ii}1/ni  (51) 

qualifying  1 
entries 

where  the  qualifying  entries  of  each  row  are  determined  from 
those  entries  where 

ABS  (C^)  >  ai  *  10“m  (52) 

where  a.  •  max  ABS  (C, .)  (largest  entry  of  ith  row) 

J  J 

m  is  chosen  by  the  user 

ni  »  number  of  C. i  entries  that  exceed  oi  *  10-ra 
threshold  J 


The  scaled  C  matrix  entries  are  then  given  by 


ILL 


ij 


Pii  QjJ 


(53) 


The  inverse  matrix  C-1  is  found  from 

Co-1  -  Q*1  C  P'1  (54) 

This  involves  evaluating  three  matrices  as  opposed  to  one 
but  the  two  diagonal  matrices  are  easily  inverted. 
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In  addition  to  the  row  and  column  scaling  techniques,  there 
are  several  perturbation  methods  proposed  in  Riference  2. 

These  involve  forming  the  matrix  A  from 

A  -  Cc  +  e  D  (55) 

where  C0  is  the  scaled  version  of  the  original  matrix  and  D  is  a 
diagonal  matrix  whose  entries  can  be  taken  as  those  of  the  diag¬ 
onal  of  A.  The  constant  e  is  chosen  to  be  suitably  small  such 
that  C  is  invertible.  (Selection  of  e  is  discussed  in  detail  in 
Reference  2.)  Then,  the  inverse  of  C  is  given  by 

C"1  =  (A  -  e  D)'1 

*  A”1  +  e  A"1  D  A"1  +  e2  (C_1D)2  Cf1  +  ...  (56) 

C-1  is  found  using  A-1,  e  and  D,  all  of  which  are  known.  A-1  was 
found  using  standard  matrix  inversion  techniques. 

_  1 

A  problem  with  this  approach  is  that  C  is  not  found  in 
closed  form  and  then  numerical  accuracy  becomes  a  question. 

The  techniques  were  exercised  using  the  computer  routine 
provided  in  Reference  2.  These  routines  were  obtained  from  Mr. 
Daniel  Kenneally  of  Rome  Air  Development  Center  who  received 
them  from  their  originator,  Dr.  V.  X.  Jain  (Reference  2). 

The  computer  routines  were  exercised  for 
response  of  a  system  with  the  linear  transfer 

_  2.8060192  x  105 

H(8)  ■  - 7f 

s  +0.011550098  (2ir)  x  10° 

_ 2.7368441  x  IQ8 

a  +  10.616086  (2n)  x  106 

The  input  to  the  system  was 

a-t  cut 

x(t)  ■  (e  +  e  )  (58) 

where  -  10?,  ao  *  1.75  x  10^  rad/sec.  The  resultant  second 
order  response,  Y2(s),  has  poles  at 


the  second-order 
function  given  by 


(57) 
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-0.011550998  (2ir)  x  10® 


'6 


'8 


'0 


'10 


'll 


'12 


-10.616986  (2ir)  x  108 
2  \t  -  -0.023101996  (2ir)  x  106 

2  X2  •  -21.233972  (2ir)  x  106 

\1  +  X2  -  -10.628537  (2ir)  x  106 

°1  +  X1  *  -1-603100429  (2 w)  x  106 

al  +  X2  “  -12.20853543  (2ir)  X  106 

a2  +  xi  -  -2.796762502  (2*)  x  106 

a2  +  x2  *  -13.4021975  (2ir)  x  106 

2ax  -  -3.183098862  (2ir)  x  106 

2a2  -  -5.570423008  (2n)  x  106 

ai  +  a2  -  -4.370760935  (2ir)  x  106 


The  C  matrix  entries  are  given  by 


C 


ij 


eXJT  i  Tn>  -  1 

17  m-1  (m  -  1)!  X1  +  1  "  m 


i  J 


(59) 

1. . . • ,12 

(00) 


T  was  set  to  600  x  10~9  second  for  this  example. 

The  computer  routines  were  exercised  by  varying  the  dimen¬ 
sion  of  the  C  matrix  from  8  to  12,  In  each  case,  the  poles  used 
were  j  »  1,..., matrix  dimension.  The  results  are  given  in 
Table  4.  The  computer  routines  evaluate  the  accuracy  of  the  in¬ 
verse  as  follows.  The  auxiliary  matrix  BQ  is  formed  as 

BQ  -  CC"1  -  I  (61) 
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where  C  is  the  original  matrix.  The  RMSE  is  defined  as 

N  N  9 

RMSE  =  E  Z  (B  y  (62) 

i=l  j*l  °ij 

which  is  the  sum  of  the  squares  of  all  the  entries  of  BQ. 

Note  that  B0  =  0  if  C~*  is  the  exact  inverse  of  C.  Also  in¬ 
cluded  in  Table  4  is  the  numerically  evaluated  determinant  of  the 
C  matrix. 


TABLE  4.  MATRIX  INVERSION  RESULTS 

C  Matrix  Root  Mean  Determinant 


Dimension 

Squared  Error 

of  C 

Matrix 

8 

0.149  x  10~7 

0.346 

X 

10-10 

9 

0.383  x  10~5 

0.124 

X 

IO-I4 

10 

0.44  x  10-3 

0.999 

X 

io-17 

11 

0.21 

0.222 

x 

10-22 

12 

1010 

0 

The  results  of  Table  4  illustrate  what  happens  as  the  matrix 
dimension  increases.  The  two  poles,  X2  and  *5>  are  nearly  equal 
which  causes  the  C  matrix  to  be  nearly  singular.  This  condition 
becomes  critical  as  the  dimension  increases  above  11.  The  results 
of  Reference  2  indicate  that  a  similar  matrix  of  dimension  12  was 
inverted.  It  appears  that  these  computer  results  were  generated 
on  a  computer  whose  dynamic  range  exceeds  the  1088  to  10"88  capa¬ 
bility  of  the  computer  used  in  this  analysis. 

Further  it  should  be  noted  from  the  work  of  Reference  4  that 
standard  matrix  inversion  techniques  were  able  to  invert  this 
matrix  when  the  dimension  was  8  or  less. 

These  techniques  offer  some  potential  relief  from  the  matrix 
inversion  problem  encountered  when  using  the  identification  tech¬ 
nique.  The  dynamic  range  of  the  computer  used  to  invert  the 
matrix,  however,  remains  the  dominant  limitation.  A  scaling  of 
individual  operations  is  perhaps  a  way  of  alleviating  this  prob¬ 
lem  but  this  is  an  extensive  process.  It  appears  at  this  point 
that  the  most  viable  way  of  increasing  the  applicability  of  the 
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identification  technique  is  to  find  methods  of  reducing  the  order 
of  the  second  order  system  response,  Y2(s).  The  primary  focus  of 
this  effort,  therefore,  has  been  to  examine  ways  of  reducing  the 
order  of  Yo(s)  without  restricting  application  of  the  identifica¬ 
tion  technique.  These  techniques  are  evaluated  in  the  following 
sections. 

C.  POLE  APPROXIMATION  APPROACH  TO  REDUCING  ORDER  OF  Y2(s) 

The  application  of  the  identification  technique  to  practical 
systems  increases  in  difficulty  as  the  ratio  of  the  highest  break 
frequency  to  the  lowest  break  frequency  of  the  linear  system  in¬ 
creases.  The  reason  for  this  is  that  the  poles  of  the  second-order 
response  include  poles  of  the  form 

+  X^  anc*  k  =  1 . N 

where  the  X^  i  =  1 ,  N  are  poles  of  the  linear  portion  of 

the  system. 

If,  for  the  system  of  interest, 

Xi  +  Xk  =  Xk 


for  some  l,  k  combination,  it  becomes  extremely  difficult  to 
solve  the  residue  equation 


R 


C-1Y 


(63) 


where  C  is  a  matrix  with  entries  of  the  form 


*JT 


'U 


i 

-  I 


xm  ♦  1<T> 


a. 


n»l  (a^) 


i  +  1  -  m 


(64) 


where  a^  =  X^  +  X^,  X^,  etc. 

If  two  of  the  a,  are  approximately  equal,  the  C  matrix  is  nearly 
singular  and  isJextremely  difficult  to  invert  using  standard 
matrix  inversion  techniques. 

To  expand  the  applicability  of  the  identification  technique, 
it  is  necessary  to  find  a  method  of  alleviating  the  computational 
problem  discussed  above.  The  approach  to  be  taken  in  this  section 
is  to  use  the  approximation 


xi  *  xj 


if  Xt  + 
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The  question  then  is  to  determine  the  impact  of  this  approximation 
on  the  identification  procedure. 


The  first  step  is  to  derive  the  functional  form  of  the  second- 
order  response,  Y2(,s),  for  this  approximation.  The  approximation 
must  be  treated  carefully  to  insure  the  correct  expression  for 
Yg(s)  is  obtained.  This  is  demonstrated  below. 

Assume  that  we  have  a  system  output  given  by 


Y(s) 


s  +  X, 


(65) 


The  corresponding  time  function  is  given  by 

-X.t  -X2t 


y(t )  =  A^e 


+  A2e 


(66) 


If  X2  *  Xf  +  e,  then 


-et  "Ai* 

y(t)  »  (Ax  +  A2  e  et)  e 


(67) 


and  for  e-MD, 


y(t)  =  (A1  +  A2)  e 


-X^ 


(68) 


Y(s)  = 


(A1  *  V 

s  +  X, 


(69) 


Therefore,  Y(s)  is  approximated  by  a  single-pole  system. 
However,  assume  that  the  system  output  is  of  the  form 


Y’(s)  = 


(S  +  X1)  (8  +  Xg) 


X2  ~  A1  +  A1  '  X2 
s  +  X^  s  +  Xg 


(70) 


The  corresponding  time  response  is  given  by 


R  -X  -  t  ~X  nt 

y * ( t )  -  t — S_  (e  -  e  *  ) 
A  2  "  *1 


(71) 


If  X2  ■  +  e,  then,  as  e-*-0 


Y  (s)  = 


(s  +  X1)< 


-X.t 

y '  ( t )  =  Bt  e  x  (73) 

In  this  case,  Y(s)  is  approximated  by  a  double-pole  system. 

The  approximation  for  X.  +  X..  =  X.  in  Y2(s)  must  be  made  be¬ 
fore  the  expression  for  Y2<s)  isJexpanded  into  partial  fractions 
The  second-order  response,  Y2  (s^  s2),  is  given  by 


Yri(S.  1  Sa  ) 


M  N  L  L 

£  ■,  1.  £1  ,  £  -  j  1  A.  . 


(Si  +  s2  -  ak  -  aR^ )(s2  -  a^^) ( Sj  -  a^)  ISJ.  +  0^ 


Expansion  and  simplification  of  Yg^^,  s2)  yields 


M  N  L  L  klk2 

Y2(s1’s2}  =  k  Z=1  k  l1  jl1  («!  +  ak  )(aj  +  ak  ) 
12  2  2 


s1  +  s9  -  2a.  "I  f 
.  1  _  Kg  _ 1 

(81  +  s2  '  0kj  -  ak2)J[(el  -  ak2) 


)(e2  -  ak  > 


(81  '  ak,)(s2  +  "  (S1  +  “i)(S2  -  ak> 


t  x  1 

(s.  +  a.  Wso  +  a, ) 


Association  of  variables  yields 


Y2(s) 


M  N  L  L  Ak  k2  1 

=  K^l  k2-l  i-1  j=l  (“7+  ak2)(aj  +  ak2>  [~ 


s  -  2a. 


S  *  S3 


s  -  2ak  s  +  (aj  -  aR  )  "  s  +  (ai  -  ak  )  s  +  (o^  +  Oj)| 
2  2  2 


(76) 

Equation  (76)  involves  products  of  transforms.  Specifically, 

A„ 


M  N  L  L 


Y0(s)  =  E  111 


klk2 


+  a.  )  (a.  +  a.  ) 


k1=l  k2=l  i=l  j=l  v“i  *k0'  “k 


s  -  (a  +  a  ) 
K1  k2 


s  -  2a. 


s  -  2a, 


(s  -  ‘“kj  +  ak2>]Cs  *  <“j  -  ak2>r  (s  (aki  f  a^)]  [s  *  (at  -  a^)] 


s  -  2a. 

_ 22 _ 

+  [s  -  <akl  +  ak  >][s  +  (ai +  (77) 

Since  ^  ak  for  any  i,  k,  and  +  Oj  f  ak  +  ak  for  any 
1  X  2 

i,j,  k^  and  k2,  Yg(s)  involves  poles  which  are  of  first  order  for 

the  Xj  +  Xj  *  Xj  situation. 

Therefore,  the  approximation  to  be  made  is 

Xi  +  X2  =  X2  for  ^2>>Xl  (78) 

in  the  expression  for  Y„(s).  Y2(s)  then  consists  of  simple-order 

poles.  The  functional  form  of  f2(s)  remains  the  same  under  the 
approximation . 

The  next  problem  to  be  addressed  is  the  generation  of  equa¬ 
tions  to  solve  for  the  unknown  A.  .  quantities.  The  general 

K1K2 

approach  is  to: 
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(1)  Identify  the  A^.,  1*1,  .  ..,  N  from  residues  of  poles 
at  s  =  -2\i 

(2)  Identify  the  Aam , 8.  f  m,  £,m  *1,  ....  N  from  residues 
of  poles  at  s  =  -  (Xj,  +  Xm) . 

(3)  Identify  the  remaining  from  residues  of  poles 

s  *  -c»i  -  1^2*  ....  N. 


The  approximation  implies  that  there  are  some  A^m,  8.  f  m,  £, 
m  «  1,  . N  that  cannot  be  identified  from  the  poles  *  £  *  ^m 
since  Xf  +  Xm  3  Xm  for  some  £,m  combination.  This  requires  the 
identification  technique  to  be  modified  to  allow  generation  of 
the  appropriate  equations  to  solve  for  the  unknown  Ak^k2 •  The 
derivation  of  a  new  technique  to  generate  equations  is  addressed 
below.  In  general,  there  are  (3N*/2)  +  (N/2)  nonzero  Ak^k2  to 
be  determined.  The  poles  of  Y2(s)  are  given  by 


s 


a^  ,  k^  —  1, 

1  2 


,N 


s 


«i  ^  a^  ,  i  1,  ...,  Li, 

2 


s  = 


-a 


i 


a 


j 


i » j  =  1, 


L. 


(79) 


First,  consider  poles  of  the  form  s  =  a^  +  ak  *  2H’ 

1  2 

£  =  1 . N.  The  terms  in  Y2(s)  corresponding  to  the  pole  at 

2X^  are  given  by  * 


Y2££(s)  “  i^1  ^  A££  (Oj  +  Xi)(ai  +  Xj)  s  (80) 

Let  the  residue  of  the  pole  at  2X.,  as  obtained  by  the  pencil-of- 
functions  method,  be  denoted  by  0^.  It  follows  that 


££ 


££ 


l  l  r 

Z  Z  , - — 

i=l  j=l  (aj 


X£)(ai 


+  X 


:>] 


-1 


£  *  1, 


N. 


(81) 


This  procedure  results  in  identification  of  N  of  the  coefficients, 
and  is  unaffected  by  the  approximation. 
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The  general  procedure  at  this  point  is  to  consider  the  poles 
of  the  form  s  =  ajc^  +  a^o  *  X£  +  Xra  where  £  f  m  and  £,m  *  1,...,N. 
Since  A^m  ■  Aj^  for  £,m  £  N,  the  terms  in  Y2(s)  corresponding  to 
the  pole  at  X^  +  Xm  are  given  by 


L  L 

,(s)  =  I  Ik. 
i*l  j=l 


°i *  “j. *  2X* 


[<aj  +  V(ai  +  Xi>(ai  +  “j  +  X£  +  V 

+  _  qi  *  a,1  *  2Xm _ 1 _ _ 1 

(ai  ~  Xm)(ai  +  V(ai  +  “j  +  X£  +  V  s  ~  H  “  XmJ 


(82) 


Let  the  residue  cf  the  pole  at  Xj,  +  *m-  as  calculated  from  the 
pencil-of-functions  method,  be  denoted  by  B 


£m' 


It  follows  that 


^£m  ^£m 


Xm> 


L  L  “  a^+a^  +  2X^ 

1=1  j=i  +  +  +  a.  *  + 


+  7 _ °i “j.  *  2X» . . _1 1 1 

ial  *  Xm><°i  +  V(“i  +  “j  +  +  Xm>J  )£  m  .  J . N 


£  f  m, £  <  m  . 


(83) 


This  procedure  results  in  identification  of  N(N  -  l)/2  of  the  co¬ 
efficients  provided  X^  +  Xm  f  Xm  for  any  m,£  = 

When  X£  +  Xm  =  Xp  is  the  situation,  specific  A£m,  £,m 
-  £  f  m,  cannot  be  identified  in  this  fashion. 

Assume  that  there  are  K  combinations  of  £,m,  £,m  =  1,...,N 
such  that  X^  +  Xm  z  X  There  remain  +  K  unknown  Ajtlk2  Quanti¬ 
ties  to  be  identified.  If  the  input  consists  of  N  decaying  expo¬ 
nentials,  the  residues  of  the  poles  s  =  -a.  -  X.,  i,j  =  1.....N 
can  be  used  to  identify  -  2K  unknown  A.  .  quantities.  Also 

k1k2 

identified  are  K  quantities  which  consist  of  sums  of  unknown 
A.  .  quantities.  This  is  illustrated  in  the  following  identifica- 
K1K2 

tion  technique  example. 
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Consider  a  nonlinear  system  whose  linear  impulse  response  is 
given  by: 


3  X ,  t 

h.  (t)  -  E  R.  e  t>0 
A  i-1  1 


(84) 


Re  { X  ±  >  <  0. 

The  problem  is  to  completely  specify  the  parameters  of  h0(t«,t0) 
where 

M  N  ak1tl  +  ak2t2 

h9(t1ft2)  =  E  E  A.  .  e  *  U(t9  -  t- ) 

212  kj-1  k2=l  klk2  1 


M  N 


ak-t2  +  ^*1 

+  E  E  A,  .  e  1  2  U(t.  -  t  ) 

k.,-1  k2'l  k1k2  x  z 


(85) 


The  parameters  are  M,  N,  a.  ,  a,  and  A,  ,  . 

K1  k2  K1  2 


The  number  of  residues,  Ak1k2 >  that  must  be  determined  is 
given  by 

W  =  N  M  -  N--^2~'  *  27.  (86) 

We  know  from  previous  analysis  that 

s  -  N  (N  -  l)2  =  3( 2 )2  =  12  (87) 

elements  of  the  set  (aki  +  ak2^  have  zero  residues.  The  elements 
of  {a^  +  ak2^  c  ^a&ik2^  whicn  have  zero  residues  have  the  form 
Xi(+  Xj  -  Xk,  i  f  k;  j  f  k.  Identifying  those  elements  of 
(akik2^  that  have  this  form  yields  the  following  zero  value 
Akik2  quantities 


A10,l  =  0 


(88) 


A82  “  0 
A92  =  0 


A93  ~  0 
A10 , 3  =  0 

We  can  group  all  the  remaining  elements  of  (a^^g)  to  specify  all 
the  distinct  elements  of  {aii.^}  .  In  this  group  Delow  we  include 
the  (ki,k2)  pairs  that  produce'the  given  sum  aici  +  a^g. 


ak  +  ak 
K1  k2 


(k1,k2) 

(4.1)  (5,2)  (6,3) 

(7.1)  (4,2)  (8,3)  (3,1)  (1,3)  (2,1)  (1,2) 

(9.1)  (10,2)  (4,3)  (2,3)  (3,2) 

(1.1) 


2X2  (2,2) 

2X3  (3,3) 

wh'.re  it  has  been  assumed  that 


(89) 


Xl  +  X2  *»  Xg 

X1  +  X3  X  3 

X2  +  *3  X2 
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We  now  demonstrate  how  the  equations  relating  to  Ak^ko  are  ob¬ 
tained  for  the  set  of  poles  at  (a^i  +  ftk2^*  Assume  that  the 
method  of  system  identification  is  applied  and  all  appropriate 
measurements  are  made.  This  means  that  a  residue  value  is  avail¬ 
able  for  each  distinct  element  in  the  set  (ajll  +  &k2^- 

We  now  consider  the  situation  that  arises  for  L  =  3.  The  result¬ 
ant  expression  for  Y2(s)  is  given  by 


10  3  3  3 

Y  (s)  =  I  E  Z  I  A.  . 

kt=l  k2=l  i=l  j-1  k1k2 


“i  +  +  2ak! 

(“j  *  ak1)(“l  *  ak1)(“l  +  +  akt  +  ak2> 


S  -  (a.  +  a.  ) 

K1  k2 


(°j  +  ak,)<ai  +  s  +  “j  - a- 


k. 


1  1 _ * 

<“i  +  akj)(aj  +  ak2^J  +  «i  -  ak2 


("i  +  °j  *  2ak2> 


+  <aJ  *  ak2)(al  +  ak2)(al  *  “j  +  akx  *  ak2>  s  +  “i  +  “j 


90) 


The  quantities  A^ko'  for  kl  =  k2 •  ^1.^2  =  are 

identified  from  the  poles  s  =  -2Xki>  ^1  =  This  identi¬ 

fies  A11(  A22,  A33. 

By  considering  each  set  of  poJes  separately,  equation  (90) 
can  be  expressed  as  three  terms 


10 

E 


3 

E 


3 

E 


3 

E  A, 


Y„  (s)  —  u  u  u  u  n,  . 

£l  kx=l  k2=l  i=l  j=l  K1K2 

ai  +  aj  +  2ak. 


(c,j  +  “k^^i +  “k^^i  *  “j  +  “kj  *  ak2>j 


s  -  (a.  +  a.  ) 

K1  „2 
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2(ax  +  a2  ♦  2ai£i) 


<*2  +  Ak  ^^ot\  +  ak  +  a2  +  Rk  +  ak, 

11  1  < 


a2  +  ak1)(2a2  +  ak 


+  ak  ) 

1  k2 


2(a1  +  a3  +  2ak^) 

(«3  »  »„  )(«!  +  »ki>(«i  +  “3  ♦  Sj  +  ak2: 


' 

2 

1^3  +  aki)(2«3  +  aki  +  ak2) 

*  m 

2  ( ot  2  ^  ®  ^  2ak  ) 

'  <03  +  akiKo2  ♦  aki>U2  ♦  o3  -  aki  * 


Y,  (5) 
Z4 


(ak  +  ak 
K1  k2 


10  3  3  3 

EE  E  E  A.  k 

ki-l  k2-l  1*1  j-1  K1K2 


“l  +  “J  +  2ak2 

<0J  *  ak2>(t,l +  ak2)(0,l  *  “j  *  ik1  ♦  ,k2) 


•  (a  +  «i  ♦  <*j) 


-  A* .  **'■»  -"***  • 
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8  ♦  a2  "  ak. 


+ °3  *  g 


:°3  *  V  (“I  * 


\  •****, 


8  *  °3  ‘  ak. 


(95) 


The  previous  analysis  has  shown  that  the  unknowns  involved  in 
equation  (95)  are 

A41’  A42’  A43*  A52 '  A63'  A71’  A83*  A91’  A10.2'  A21’  A31'  A32 

If  we  consider  the  pole  at  s  »  -  and  follow  the  pro¬ 

cedure  detailed  previously,  we  obtain  the  equation  from  (95) 

A  1  /  1  +  1  ,  1  \ 

"41  a,  +  a.  la,  +  a,  a _  +  a,  a „  +  a.  I! 

al  a4  \al  *1  a2  al  a3  al/J 

+  4  _ l _  /  1  +  1  +  1  \1 

h21  aj  +  a2  lc^  +  a2  +  al  a3  +  al/ 

+  Anj^  -!  ♦  a7  (“T^l  +  a2  ♦'  *1  *  “3  *  ‘i)] 

*  *3lf  SJ  i  .3  («!  *  *J  *  Bj  *  *1  +  “3  *  al)j 

+  A  [ _ l _  / - 1 -  +  - i -  +  i \|  -  D,  (96) 

*91  Bj  ♦  -9  (bj  ♦  a,  »2  *  »J  b3  *  a^j  1 
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i 


which  reduces  to 


A41  +  A21  A71  A31  +  A91  „  5 

°1  +  a4  al  +  a2  “l  +  a7  al  +  a3  al  +  a9  1 


(97) 


Similarly  the  pole  at  s  ■  “  a2  +  ai  yields  the  equation 

A, 


r  i 

l  *  A  r 

1 

i . .  r  i  i 

L<°2  +  a4> 

J  21  L 

°2  +  a2j  “71La2  +  a7J 

r  i  i 

+  A  r 

1 

JL  ft  • 

11  L°2  +  a3J 

A91 

>a2  +  a9^j  2 

>  pole  at  s 

"  -Q3  +  al 

yields 

the  equation 

f  1  _  1 

+  *  r 

1  1 

+  A  r  1  1 

l<°3  *  a4>J 

21  l°3 

+  a3J 

71  [<a3  +  a7>J 

r  1  l 

*  A  1  — 

1  .1 

m  R  1 

!1n  *  “a) 

01  L“3 

+  a9  J 

D3 

(98) 


(99) 


From  the  previous  analysis 


a4  -  0 


a2  “  X2 


a7  "  Xg  -  XX  a^  ■  Xg 


a9  X3  "  X1 


(100) 


From  the  constraints  on  the  analysis,  it  is  known  that 
a4  *  a7  +  a0 

and  that  our  approximation  requires 
a9  "  a3’  a7  *  a2 
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Equation  (100)  reduces  equations  (97)  and  (99)  in  matrix  form,  to 


■ 

9  m 

»  - 

1  1 

“l  +  a7  “l  +  a9 

A41 

V 

1  1 

— -  t 

“2  +  a7  °2  +  a9 

A71  +  A21 

S3 

D2 

1  1 
a3  +  a7  a3  +  a9 

m 

A91  +  A31 

- 1 

CO 

IQ 

_ i 

(101) 


Previous  analysis  (Reference  3)  has  shown  that  this  set  of 
equations  is  linearly  independent  and  therefore  can  be  solved  for 


A41’  ^A71  +  A21^  and  (A91  +  A31^ 


If  this  process  is  continued  by  considering  the  pole  at  s  ■  -a- 
+  a o,  an  equation  similar  to  equation  (97)  is  obtained,  namely1 
with  the  a-  quantity  replaced  by  a-2.  Similarly,  this  occurs  for 
the  equations  (98)  and  (99).  Under  these  conditions,  equation 
(101)  becomes 


‘42 


A52  +  A10 , 2 


12 


‘32 


V 


D6 ' 


(102) 


If  L  ■  4  in  this  case  the  above  set  of  equations  would  be  linear 
ly  independent  and  solvable  for 


A42 '  A12 ’  A32 ’  and  ^ A52  *  A10,2^ 

Knowing  A1<3  (also  A01 )  permits  identification  of  A71  from  previous 
analysis. 
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By  considering  the  pole  locations  s  ■  -a^  +  ag,  i  ■  1,  2,  3. 
the  following  matrix  equation,  similar  to  (10i)  and  (102)  is  ob¬ 
tained  : 


’ll  1  1 

A 

m  m 

n  1 

«1  al  +  a6  al  +  a8  “l  +  a3 

*43 

°7 

11  1  1 

A 

a2  a2  +  a6  “2  +  a8  a2  +  a3 

*63 

S 

D8 

11  1  1 

A 

_  1 

n 

°3  a3  +  a6  a3  +  a8  a3  +  a3 

*83 

L  9  J 

A13 

If  L  -  4  in  this  case,  these  equations  are  linearly  independent 
and  solvable  for 

A43’  A63’  A83 ’  A13 

A  problem  may  arise  if  »i  +  a3  =  at  -  ag  ( ag  *  -a3>  for  all  i. 

If  this  does  not  occur,  then  all  the  Aj^kg  quantities  have  been 
identified  except  for  A52,  A^q  «.  The  sum  A52  +  Aio?2 
has  been  identified.  However, ’  A52  can  be  found  from’the  residue 
of  the  pole  at  s  *  -Xj..  This  portion  of  the  response  is  given  by 
equation  (91).  The  only  unknown  in  this  equation  is  As« • 

Once  A52  is  known,  Aiq,2  is  also  known.  Therefore  all  the  unknown 
Akik2  Quantities  have  6een  identified. 

The  problem  at  this  point  is  to  demonstrate  that  this  can  be 
achieved  for  all  N. 


We  now  attempt  to  develop  a  technique  which  works  in  general 
for  all  values  of  N.  For  convenience  we  order  the  X*, 
i  -  l,...,N,of  the  linear  impulse  response  such  that 

*1<*2<*3<*4,’*<XN 


The  identification  technique  will  yield  a  set  of  equations  at  pole 
s  s  _0j  +  Xj  given  by 


N+K' 

t 

s=l 


1 

a 


k 


Is 


(104) 
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where  Di  is  the  residue  of  pole  s  ■  -oj  +  Xj  with  all  known 
quantities  subtracted  out.  The  index  R^g  corresponds  to  those 
values  of  kj  for  which  is  unknown. 

The  N  +  K'  akjs  Quantities  are  Riven  by 


XN  “  X1 

=  Xj  K’  of  these  values, 

j  =  N,  N  -  1 . N  -  K’  +  1. 


The  entries  correspond  to  those  poles  that  arise  because 
of  the  pole  approximation.  For  example,  if  X„  +  X,  =  X.,  then 
X-  *  a,  in  the  above  set. 
is 

Similarly  for  a 2,  we  obtain 


kls2 


N+K" 

s-1  <iTTs“)  ~  ~2 


where 


K"  <  K' 


and 


(106) 


X„  -  X„ 
N  2 


Xj  )  K"  of  these  values,  J  *  N,  N  -  1 . N  -  K"  -  1 . 


\ 
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Continuing  this  process  yields 


N+K1* 

s-1  (ctj  +  ak 


klsN 


N 


Is 


where 

KN  <  KN_1  <  ...  <  K"  <  K' 

and 


a 


k 


Is 


x 


2  " 


X 


N 


(107) 


XN-1  “  XN 

N-1,...,N-KN+1 

Solution  of  these  equations  will  result  in  identification  of 
a  number  of  the  Aj^^o  quantities.  However  several  coeffi¬ 

cients  will  combine  together,  and  only  a  sum  of  coefficients  is 
obtained.  This  is  illustrated  below. 


Assume  that  K'  =  1.  This  implies  that 

*1  +  ~  XN  and  X1  +  *i  Xi  for  1  =  1»  ....  N  -  1. 

The  equations  derived  are 


N+l 

E 

s*l 


Qt  . 

J 


"is1 


+  a, 


‘Is 


(108) 
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where 


X 


1 


XN-1  “  X1 
XN  "  X1  = 


N 


If  L  *  N,  then  a  set  of  linear  independent  equations  is  gener¬ 
ated  that  results  in  solution  for  all  the  A.  -  quantities  except 

for  a  pair  given  by 


For  k2  *  2, 


N  -N+3, 1 
A 


+  A 


N1 


kls2 


N 

s«l  aj  +  ak 


Is 


(109) 


where 


kls 


*  X1  -  X, 


XN  "  X2 


These  equations  are  linearly  independent  and  can  be  solved 
directly  for  A,,  0.  Similarly  a  solution  is  obtained  for  the 

Ki 

quantities  A.  A,  w  -  . 

Is*3  klsw”1 


For  kg  *  N,  the  equations  are 
A, 


klsN 


N+l 
s*l  aj  +  ak 


=  D 


N 


(110) 


Is 


52 


where 


XN-1  '  X1  “  XN-1 
’  XN  -  X1  *  XN 

Solution  of  these  equations  leads  to  identification  of  N  -  4 
A.  .  quantities  and  the  summations 
k1k2 
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AN-1,1  +  AN2-N+2,1 

AN,1  +  AN2-N+3,1 
For  k  =  2,  the  equations  are 

?  Ss2  ,D 

s-1  <“j  +  ak  >  2 


(112) 


These  are  solved  directly  for  the  Ak^2'  Similarly  for 
k2  -  3,...,  N.  This  implies  that  A1|N  and  Ai  n-1  are  identified 
provided  N  >  2.  Using  the  symmetry,  Afcjj^  *  «k2ki>  kl>k2“  1*  •••»”> 
yields  solution  for  all  the  Akik2- 

Assume  that  K'  *  2  but,  in  this  case,  the  poles  are  such  that 
\N  +  *i  a  XN  and  XN  +  X2  -  XN<  The  resultant  equations  are 


N:2  ^is1 

s-1  “j  +  ak. 


(113) 


where 


X2  ~  X1 


XN  ~  X1  =  XN 


tion 


This  yields  solution  of  N  -  2  Akj>sl  quantities  and  the  summa- 


Vl  +  AN2  -N+3,1 


For  k2  *  2,  the  equations  become 

N+2  Ak.  2 
E  .  -  -  ls - ,  *  D9 

-1  (“j  °kis 


(114) 
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where 


kls  *  0 


xi  "  x2 


X 3  ■  X2 


XN  *  X2  =  XN 

Solution  of  these  equations  yields  all  of  the  unknown  Akjk2 
except  for  the  summation 

AN2  +  AN2-N+4 , 2 


For  k2  *  N,  the  equation  becomes 

A,. 


N-l 

l 


_ =  D 

s-1  aJ  +  Akls  N 


(115) 


where 


*  0 


k  Is 


X1  "  XN  '  "XN 
X2  -  XN  2  "XN 


XN- 1  "  XN 


The  solution  of  these  equations  identifies  all  Ak  quanti¬ 
ties  except  for  the  summation,  KlsN 


A2N,N  +  A3N-1,N 

This  yields  solution  for  A^  which  implies  that  the  only  unknown 
Ak1k2  at  this  point  are 

A2N,N  and  A3N-1,N 
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An  equation  involving  A2N  and  As^-i.N  can  obtained  from 
the  residue  of  the  pole  at  s  *  -2ai.  This  equation  has  the  form 


(«t  »  *  (ai  ♦  ‘  “i  <116> 

where  *  residue  at  pole  at  s  ■  -2a^.  This  equation  can  be 
rewritten  as 


3N-1.N 


A2N , N  .  A3N-1 , N  u  , 
Jo1  +  X^  2a  i  +  Xg  1 


(117) 


The  two  equations  involving  AgN  N  and  N  can  be  repre¬ 

sented  in  matrix  form  as 


k2N ,  N 


(118) 


2a1+  X1 


2al  +  ^2 


l3N-l ,  N 


[B]  [A]  «  [Z] 

For  linear  independence 
det  B  f  0 

For  det  B  =  0,  it  is  required  that 


(119) 


(120) 


2a^  +  Xg  2a^  + 


(121) 


or  X0  =  L,  which  is  not  permitted  by  assumption. 

a  1 

Therefore,  it  is  possible  to  solve  for  all  the  A^  ^  quantities 
even  if  K‘*  Z.  *  ^ 


Further  analysis  has  not  led  to  a  method  of  demonstrating 
that  the  Aklk9  coefficients  can  be  identified  for  an  arbitrary 
number  of  pole  pair  approximations  (X1  +  Xj  =  X  ^ ) .  One  of  the 
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reasons  for  this  is  illustrated  in  the  analysis  for  two  pole  pair 
approximations  (K'  =  Z).  There  are  many  combinations  of  potential 
pole  paires  to  achieve  K'  =  2.  As  K1  increases,  these  combina¬ 
tions  increase  significantly  in  number,  and  it  is  necessary  to 
show  that  a  set  of  linearly  independent  equations  can  be  found. 
This  process  rapidly  becomes  complex  to  show,  in  general ,  that 
these  equations  can  be  generated. 

At  this  point,  it  should  be  noted  that  significantly  more 
equations  are  generated  than  are  used  in  the  identification  pro¬ 
cess.  It  is  possible  that  some  of  these  equations  can  be  used  to 
form  a  linearly  independent  set  of  equations  to  solve  for  the 
unknown  Aj^^,,-  The  problem  is  that  a  method  has  not  been  found 
for  der onstrating  that  a  linear  independent  set  of  equations  in¬ 
volving  the  unknown  A^iko  quantities  can  always  be  obtained  using 
the  identification  technique.  The  algebraic  nature  of  the  equa¬ 
tions  has,  to  date,  prevented  linear  independence  of  any  set  of 
equations  other  than  those  generated  from  poles  at  8  *  -a*  +  Xk2> 
i  «  1 . L;  k2  *  1,...,N  from  being  proved  in  general. 

The  above  analysis  implies  that  the  pole  approximation 
(Xi  +  *  Xj)  is  a  reasonable  approach  to  alleviating  the  nu¬ 

merical  problem  associated  with  poles  of  this  type.  For  one-  or 
two-pole  pairs,  it  has  been  shown  that  the  Aj^j^  quantities  can 
be  identified.  For  more  than  two  poles,  each  system  identifica¬ 
tion  problem  must  be  considered  individually.  Although  a  general 
identification  technique  is  not  defined  here,  it  is  probable  that 
there  will  be  a  sufficient  number  of  linear  independent  equations 
to  solve  for  all  the  Afcik2  quantities.  Each  system  identifica¬ 
tion  problem  must  1 2  considered  individually  to  find  a  sufficient 
number  of  linearly  independent  equations. 

D.  DOMINANT  POLE  CONCEPT 


As  noted  previously,  the  identification  technique  becomes 
significantly  computationally  complex  as  the  number  of  poles  of 
the  linear  system  increases.  The  second-order  impulse  response 
of  a  weakly  nonlinear  system  is  given  by 


^2^1  ’  ^2  ^ 


kl=1  V1 


Ak  k  e 
k1k2 


Sjh  *  **2X2 


U(t2  -  tx) 


ak  tl  +  ak  t2 
2  1 


M 

+  E 


N 

E  A,. 


e 


U(t,  -  t0)  (122) 


For  a  given  system,  it  is  possible  that  this  impulse  response 
is  dominated  by  a  limited  number  of  terms.  This  implies 

that  h2<t^,to)  can  be  represented  by  fewer  terms  than  are  given 
in  equation  (122).  This  reduces  the  number  of  coefficients  that 
must  be  identified  and  will  ease  the  computational  problem 
associated  with  the  identification  technique.  This  section  in¬ 
vestigates  how  thi9  approach  may  be  implemented  to  reduce  the 
order  of  the  second-order  impulse  response. 

The  second-order  system  response,  Y2(s),  has  been  shown  to 
be  given  by  equation  (19).  The  primary  focus  on  reducing  the 
order  of  the  identification  problem  must  be  on  reducing  the  order 
of  Y2(s).  This  is  because  the  order  of  the  residue  equation 
R  -  C-l  Y  is  directly  determined  by  the  order  of  Y2(s). 

All  of  the  poles  of  Y2(s)  do  not  contribute  equally  to  the 
system  output,  y2(t).  It  is  possible  that  some  poles  of  Y2(s) 
contribute  negligibly  to  the  output.  If  these  poles  can  be 
identified,  then  the  second-order  lesponse,  Y2(s),  can  be 
approximated  by  a  response  Y^(  s) ,  where  Y2(s)  does  not  contain  the 
poles  of  Y2(s)  that  negligibly  impact  the  output. 

The  basic  problem  with  identifying  which  poles  of  Y2(s)  have 
a  negligible  impact  on  the  output  is  that  this  requires  knowl¬ 
edge  of  the  Akik2  quantities  that  we  are  trying  to  identify. 

This  implies  that  the  poles  of  Y2(s)  that  have  a  minor  impact  on 
the  system  output  must  be  identified  by  other  than  analytical 
means.  However,  it  is  unlikely  that  this  can  be  accomplished. 

A  potential  method  of  reducing  the  order  of  Y2(s)  by  identi¬ 
fying  the  negligible  poles  is  presented  here.  The  number  of 
poles  and  the  location  of  each  pole  of  Y2(s)  are  known  from 
identification  of  the  linear  transfer  function,  Hi(s).  The  pro¬ 
posed  technique  is  as  follows.  First,  the  nonlinear  system  is 
excited  by  the  input  x(t)  »  e"®!1.  The  resultant  output,  Y2(s), 
contains  N[(N  +  5)/23  +  1  poles.  The  normal  procedure  is  to 
integrate  the  input  and  output  N'  times  and  solve  the  result¬ 
ant  equations  for  the  system  residues.  The  suggested  procedure 
is  to  assume  that  the  number  of  significant  poles  of  Y2(s)  is 
N"  <  N'.  The  input  and  output  responses  are  integrated  N"  times, 
the  appropriate  inner  products  are  formed  and  the  resultant  equa¬ 
tion  for  the  pole  locations  is  solved.  If  these  pole  locations 
agree  with  N"  of  the  predicted  pole  locations,  it  can  be  assumed 
that  an  N"  pole  approximation  of  Y2(s)  is  a  valid  representation. 
If  the  pole  locations  do  not  agree,  then  N"  is  increased  by  1 
and  the  above  process  is  repeated.  This  is  continued  until  there 
is  good  correlation  between  the  identified  poles  and  those  pre¬ 
dicted  from  H^s). 
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Id  order  to  demonstrate  that  such  &  technique  is  feasible, 
consider  the  following  example.  Suppose  we  have  a  system  with 
three  poles,  where 

H(s)  »  2.8069192  x  105 

s  +0.011550998  (2tr)(106) 

2.7368441  xlO8 


s  -  10.6161986(2tr)(106) 

r3 

+  - 2 -  (123) 

s  -  1.25(2ir)(10°) 


The  residue,  83,  will  be  left  undefined  for  the  moment.  This 
transfer  function  was  inserted  into  the  computer  simulation  of 
the  identification  technique.  (This  simulation  was  described  in 
detail  in  Reference  1.) 

The  simulation  assumed  that  the  system  of  interest  was  a  two- 
pole  nystem.  The  residue  value,  R3,  was  varied  in  amplitude  to 
determine  under  what  conditions  H(s)  could  be  accurately  repre¬ 
sented  by  a  two-pole  system. 

The  results  of  the  simulation  are  given  in  Table  5.  These 
results  indicate  that  for  R3  <.  104 ,  identification  technique 
identifies  the  other  poles  and  residues  of  H(s)  with  less  than 
0.25  percent  error.  For  these  cases,  R3/2. 8069192  x  105  <  0.035 
and  R3/2. 7368441  x  108  <  3,6  x  10"5,  which  indicates  that  the 
poles  of  H(s)  at  s  =  -0  11550998 (2tt)  x  10®  and 

s  *  -10 . 616986( 2tt )  x  10b  dominate  the  transfer  function.  As  the 
residue  Rq  approaches  and  exceeds  the  magnitude  of  the  residue 
of  the  pole  at  3  =  -0 . 01150998( 2ir ) ( 106)  ,  the  identification  tech¬ 
nique  degrades  in  performance  as  it  attempts  to  identify  the  two- 
pole  model  of  H(s).  The  reason  for  this  is  that,  as  R3  increases, 
H(s)  is  not  dominated  by  the  two  poles;  therefore,  the  proposed 
technique  is  no  longer  a  valid  approach  to  identifying  H(s). 


Another  approach  to  reducing  the  order  of  Y2(s)  and  subse¬ 
quently  easing  the  computational  problem  is  to  apply  a  dominant 
pole  concept  to  the  linear  portion  of  the  system.  The  linear 
transfer  function  of  the  system  of  interest  is  assumed  to  be  of 
the  form 


H1(s) 


It  is  possible  that  Hi(s)  is  dominated  by  several  of  the  N  poles. 
For  example,  suppose 
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TABLE  5.  IDENTIFICATION  TECHNIQUE  PERFORMANCE  -  TWO-POLE  APPROXIMATION 

OF  A  THREE  POLE  SYSTEM 
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H1(s) 


s  -  Xj  s  •  \ 2 

r  _  a_vu^2_\ 

8  '  R  +  R  ’ 

"  (R1  +  *V  (S  -  Xx)  <8  -  \2) 


(124) 


or  if 


R1  X2  +  R2  X1 

Rx  “  R2  "  X1 


R1  X2  +  R2  X1  \n 
R1  +  R2 


(125) 


then  H1(8)  could  be  represented  by 


lys)  = 


(Ri  ♦  r2> 

8  -  X„ 


Hx(s) 


(R1  *  V 

8  -  A, 


(126) 


respectively.  In  these  cases,  the  two-pole  linear  transfer  func¬ 
tion  Hi(s)  can  be  represented  by  a  single-pole  transfer  function. 

This  is  significant  because  the  order  of  Y2(s)  varies  approxi¬ 
mately  with  N^,  so  that  any  reduction  in  N  achieves  an  even 
greater  reduction  in  the  number  of  poles  of  Yg^). 

The  impact  of  this  approximation  on  H2(si,s2)  or  h2(ti,t2) 
is  presented  below.  Consider  an  exact  system  representation  given 
by 


Hi<8>  "  1  rrr 

x  i-i  8  Ai 


(127) 


iikLiji 


The  resultant  h2(ti,t2>  is  given  by 


x1(t1  +  tg) 

’  A11  e 


+  A12  e 


X1  tl  +  X2  *2 


a  oX2  *1  *  Xl  >  »  X2^1  +  t2) 

"21  e  +  a22  e 


+  A  e  1  "2  ♦  A  e*2*2 
"31  "32  6 


+  A42  e 


+  A51  e 


X-  t, 

1  4 

(Xx  ~  Xg)ti  +  Xg  tg 

(X5  —  X^)t^  ^  X ^  tg 


U(tg  -  tj) 


X^(tj  +  tg)  X^  tg  +  Xg  t 

An  6  +  A12  6 

X2  t2  *  X1  tl  .  X2^tl  +  t2 
+  A21e  +  A22  e 


If 


then 


, .  X1  tl  ,  .  X2  tl 
+  A31 e  +  A32  e 


+  A42  e 


(xx  -  x2)t2  +  x2  t1 


*  A61  e 


( X 2  ~  ^j)^2  *  ij. 


U(tx  -  t 


R1  X2  +  R2  X1  >  , 
Ri  +  fto - X2 


Hi<8>  “  «nrr: 


The  resultant  hg(t1(tg)  is  given  by 
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1 

) 


2)  (128) 


r 


•,  A'  A*N.  ‘f  fi, 


+  t2> 


^2^1  ’  ^2  ^ 


+ 


U(t2 

U(tj 


V 

t2)  (129) 


The  approximate  expression  for  hgCt^tg)  implies  that 

A11  '  A11 

A21  *  A31 

A12  *  A21  =  A22  =  A32  =  A42  =  A51  =  0  (130) 

for  h2(t1,t2)  ■  h2(t1,t2). 

The  approximate  expression  for  h2(ti,t2)  retains  only  two  of 
the  original  Akikg  coefficients,  meaning  that  the  approximation 
(Rl  ^2  +  r2  ^l)/(Rl  +  r2)  —  *2  implies  that  six  of  the  Akik2  co¬ 
efficients  are  zero.  In  order  to  see  how  this  approximation 
comes  about,  we  consider  the  following  example,  which  is  a 
simple  single  nonlinearity  (no-rnemory)  nonlinear  system,  as 
shown  in  Figure  4. 


a 

L  K„  v" 

n»1  n 


Figure  4.  Simple  Single  Nonlinearity  (no-memory)  Nonlinear  System 

It  has  been  shown  (Reference  3)  that  the  second-order  impulse 
response  of  this  system  has  the  form 
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*1  *  ^2> 


N  N 
E  E 

i“l  J“1 


/  RlRiRj 
\^1  ‘  Xi  “ 


~  X  ^ )  t  ^  + 


X1  *2 


-  e 


Vi + 


l2  <  *1 


where  it  has  been  assumed  that  the  linear  impulse  response 
system  is 

N  X  t 

h1  ( t )  =  E  R.  e  1  t  >  0  . 

i  =  l  1 

For  purposes  of  example,  it  is  assumed  that  N  *  2  and  that 

X.  t  X0  t 

hx(t)  -  R1  e  1  +  R2  e  z  . 

If  h2(ti,t2)  is  expanded  for  N  =  2  and  put  in  the  standard 
tional  form  given  by 


5  2 


ak  tl  +  ak  t2 

h9(t1(t9)  =  E  E  A.  k  e  1  2  U(t0  -  t. ) 

Z  1  *  k1=l  k2=l  klk2  2  1 


ak  +  f  2 

k  k  e  2  1  U(tl  -  *2 

k1=l  k2=l  K1K2  1  4 


5  2 

+  E  E  A 


where  the  A.  .  quantities  are  defined  as  follows: 
K1K  2 


(131) 
of  this 

(132) 

(133) 
f  unc- 


) 

(134) 
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ti 

F  ! 


lll 


K  R  2  r  h  +  R2 

21  L  X1  X2  -  2Xd 


A12  "  A21 


■-  w2[^y 

■  w2  [-xf  +  r^r2\ 

2  r?i  t  h  1 

1  LX1  X2  J 

[!}-5] 


22 


A31  =  K2R 


A32  "  K2R2 


A  41  =  0 


A42  ”  K2R 


‘51 


A  [x,  -2X,  *  x‘] 
-  K2Ria2  [>;v-2>1  +  r[] 


A52  =  0  (135) 

The  dominant  pole  assumption  was  that 
R-.  X9  +  R9  X- 

Rx  +  R2  X 2  '  (136) 

This  requires  that 

R1  >>  R2  an<*  R1  X2  >:>  R2  X1 

If  this  assumption  is  applied  to  the  resultant  quantities, 

we  obtain 


!L 


X 

k 

'  t 


which  reduces  the  second-order  impulse  response  to: 

r  kx3  h  +  Xj  t2 

h*(ti-*»>  '  K2  —  e 


«(t2  -  V 


Ri3  xi  *»1 

+  "e  J 

r  Rj3  xx  t2  ♦  xx  tj 
♦  K*  [-  — 6 


bi  xi  ii 

*  IT  * 

X1 


U(t2  -  tx) 
If  h^t)  is  approximated  by 


(139) 


h1(t)  =^6 


Xe  t 

1  t  >  o 


(140) 


uslng  the  dominant  pole  concept,  the  resultant  «*«,.»,>  i=  S-n 


by 


,  K2  El3  feXl  -  eH 

h2(t1.t2>  ‘  l® 


or 


[' 


+  'Au  e 


X1  *2  +  xl  *1  + 


1  +  X2  l2  ^ 

U(t2  -  tj) 

*2  +  X2  4lj 

o(tj  -  V, 

X1  l2 

A21  8 

]  VI* 2  '  tl) 

X1  t] 

+  A21  6 

LJ  D(tt  -  t2 

(141) 


(142) 


where 

A 


11 


K2  V 


*2  V 


(143) 
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(144) 


It  is  noted  that,  given  the  dominant  pole  assumption, 

This  analysis  demonstrates  how  the  approximation  in  hi(t) 
propagates  to  h2(ti,t2).  It  supports  the  approach  using  the 
dominant  pole  approximation  on  the  linear  transfer  function. 

It  is  noted  from  the  analysis  that  the  dominant  pole  is  pres¬ 
ent  because  H^(s)  has  a  pole  and  zero  which  tend  to  cancel  each 
other  in  the  s  =  o  +  Jco  domain.  This  illustrates  why  this  approach 
cannot  be  used  on  Y2(s)  since  the  unknown  Au.^2  quantities  pre¬ 
vent  the  zeros  of  Y2(s)  from  being  known.  Tnis  prohibits  associ¬ 
ation  of  pole-zero  pairs  for  possible  cancellation  and  subse¬ 
quent  reduction  in  the  order  of  Y2(s). 

E.  RESTRICTED  FREQUENCY  RANGE  CONCEPTS 

There  are  several  methods  cf  reducing  the  order  of  the  second- 
order  response  which  we  o  irissif  at-  restricted  frequency  range 
approaches.  These  approaches  busically  modify  the  input  or  the 
system  output  t  o  ease  the  identity oat  ion  problem. 

The  primary  restricted  frequency  approach  is  to  use  a  filter 
on  the  output  of  the  system  under  tes': ,  The  purpose  of  the  fil¬ 
ter  is  tc  selectively  restrict  f he  sy.uem  output  to  a  particular 
frequency  ranee.  Consider  the  example  shown  below: 


NONLINEAR**!  j  p, 

1  SYSTEM  J - *  1 


(0 


filter 


For  second-order  impulse  response  considerations,  the  equivalent 
system  is  as  shown  below: 


V2(s)  F,(s) 
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The  system  output  becomes  Y2(s)  F^s).  Fi(s)  must  be  band- 
limited  with  respect  to  Y2(s).  The  number  of  poles  of  Y2(s)Fi(s) 
depends  on  the  number  of  poles  of  Yo(s)  within  the  bandwidth  of 
F-(s)  and  the  number  of  poles  of  F^(s).  For  this  technique  to 
offer  any  advantage,  the  number  of1 poles  of  Y2(s)  F1(s)  must  be 
less  than  were  present  before  filtering. 

The  poles  of  Y2(s)  are  of  the  general  form 


xi 

i  = 

1, 

.  .  N 

x* 

+ 

X 

m 

2. ,  m 

= 

1, 

.  .  .  ,N 

-°i 

+ 

X  . 

J 

i  ,  j 

1, 

.  .  .  ,N 

-a  . 

1 

+ 

i.  j 

= 

1, 

.  .  .  ,N 

The  ai  are  selected  by  the  identification  technique  user.  Proper 
selection  of  the  cxi  can  cause  the  poles  of  Y2(s)  to  bunch  up  in 
certain  frequency  ranges. 


Consider,  for  example,  a  two-pole  system  with  poles  and  X2 
(X2  >  A!).  The  poles  of  Y2(s)  are 


X 


1 


-2a 


1 


X 


2 


-(«i  *  a2) 


A  +  A  -2a 

12  2 

-a  +  X  2A 

11  1 


-al  +  X2 


2X 


2 


-a 


2 


+ 


X 


1 


-a 


2 


+ 


X 


2 


Suppose  that 


(145) 
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Then  the  poles  of  Yg(s)  are 
X1 

P1X1 

2\x 

2p1x1 

(1  +  p1)x1 
(P2  +  1)X1 

(p2  +  P1)X1 
(p3+  i)x1 

(  P3  +  P1^X1 

2p2xi 

(p2  +  p3)x1 

2p3Xx 

Assume  that  we  select  p.  *  4.5,  p2  =  1.2,  p„  =  3.9,  then  the 
poles  are  J 

\x,  2Xt,  2 . 2X  x ,  2.4X1,  4 . 5X  x ,  4.9X1,  5.1X1(  S.SX^  5.7X1; 
7.8XX,  8.4X1,  9\1. 

These  poles  are  bunched  in  essentially  two  groups: 

I 

(Xj^  -*•  2.4X1)  and  (4.5XX  -*•  9X1). 

If  a  filter  with  a  3-dB  bandwidth  of  approximately  2.4X-  is  used 
on  the  output  of  the  system  under  test,  the  resultant  output 
contains  the  contributions  of  only  a  limited  number  (in  this  case, 
four)  of  the  poles  of  Y2(s). 
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The  performance  of  this  approach  was  evaluated  using  the 
computer  simulation  of  the  identification  technique.  The  system 
considered  was,  once  again,  that  given  by 


H(s) 


2.8069192  x  10' 


2.7368441  x  10 


8 


.6 


s  +  0 . 011550998 ( 2tt ) ( 10  )  s  +  10. 616986 (2tt  )  x  10" 

(146) 

The  filter  transfer  function  was  assumed  to  be 

7 

«A(S)  -  -  (147) 


S  +  Yi 


where  the  pole  location  y-  was  left  variable.  The  simulation  was 
to  evaluate  the  poles  of  a  two-pole  system  representation  of 
H(s)Ha(s),  where  the  two  poles  to  be  identified  are  the  two  low 
frequency  poles,  -0. 011550998(2 *)  x  10®  and  y, .  The  pole  of  in¬ 
terest  is  that  at  s  E  -0.011550998(2*)  x  10$  since  yj  will  be  * 
selected  by  the  user  and  will  be  known. 


input 


This  procedure  was  simulated  for  y^  =  0.51  (2ir)  x  106  and  an 


x(t) 


-0. 2( 2tt )  x  106t 


The  results  are  tabulated  in  Table  6.  These  results  indicate  that 
this  approach  produces  acceptable  performance  if  the  integration 
time  is  increased  above  that  used  for  the  original  system  identifi¬ 
cation.  This  i'i  to  be  expected  because,  in  this  case,  the  identi¬ 
fication  technique  is  attempting  to  identify  two  low  frequency 
poles  instead  of  one  low  and  one  high  frequency  pole;  this,  in 
general,  will  require  longer  integration  times. 

This  procedure  was  repeated  for  an  input  given  by 
x(t)  -  e'2’  x  1<,6t 

The  results  are  presented  in  Table  7.  These  results  basically 
agree  with  those  of  Table  6  and  support  the  need  for  increased 
integration  time. 
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The  identifier  tea  technique  can  then  be  used  to  evaluate  the 
residues  of  the«e  poles.  Once  these  residues  are  known,  the 
technique  can  be  repeated  without  using  the  filter.  The  known 
portion  of  Y2(s)  can  be  subtracted  out  before  processing  and  the 
resultant  identification  problem  is  reduced  to  one  of  lesser  order 
(in  this  case,  8  instead  of  12). 

An  alternative  approach  at  this  point,  once  4  of  the  12 
residues  have  been  identified,  is  to  attempt  to  select  the  to 
separate  the  remaining  poles  into  distinct  groups  and  use  a  dif¬ 
ferent  filter  to  limit  the  number  of  poles  of  Y-(s).  This  essen¬ 
tially  repeats  the  original  identification  technique  approach 
but  on  the  reduced  order  system. 

The  key  to  the  technique  is  to  use  a  filter  which  effectively 
attenuates  the  contributions  of  the  poles  outside  the  frequency 
band  of  interest.  Consider  the  example  shown  below.  The  system 
response  is  as  shown  in  Figure  5. 


Figure  5.  Example  System  Frequency  Response 


The  filter  should  have  a  break  frequency  equal  to  or  slightly 
greater  than  fL-  The  amplifier  response  should  be  down  con¬ 
siderably  at  frequency  fjj  to  effectively  attenuate  the  frequency 
response  of  Yo(s)  above  f h •  An  attenuation  of  at  least  20  dB 
seems  reasonable  for  adequate  performance  of  the  technique. 

Another  approach  to  reducing  the  order  of  Y2(s)  for  identifi¬ 
cation  purposes  is  presented  here.  It  has  been  shown  (Reference 
2)  that  the  pencil-of-functions  identification  technique  can  be 
modified  to  divide  the  frequency  spread  of  the  system  output 
into  three  bands:  low,  midrange  and  high.  The  identification 
technique  is  applied  by  selecting  an  appropriate  input  frequency 
fairly  well  matched  to  one  of  the  frequency  ranges  given  above. 
This  procedure  is  repeated  for  each  frequency  range.  The  total 
transfer  function  is  obtained  by  matching  the  functions  at  the 
transition  points  between  frequency  ranges  and  slightly  modifying 
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pole  locations  and  gain  constants.  This  technique  has  been  shown 
(Reference  2)  to  produce  accurate  results  while  reducing  the 
order  of  the  identification  problem. 


Appropriate  selection  of  integration  time  for  the  identifica¬ 
tion  technique  can  be  used  in  special  cases  to  reduce  the  order 
of  the  identification  technique.  Consider  a  two-pole  system  with 
poles  Xi  and  X2,  where  X2  >>  Xj_.  If  a  short  integration  time  is 
used,  e.g.,  T  =<  l/\\,  the  system  output  will  be  impacted  very 
little  by  the  low  frequency  pole.  This  is  because  the  contribu¬ 
tion  from  the  pole  s  =  \1  is  essentially  constant  over  the  inte¬ 
gration  period.  This  implies  that  the  variation  of  the  system 
response  over  the  integration  period  is  due  to  the  high  frequency 
pole  only. 

a  technique  would  perform, 

_ 2,7368441  x  108 

s  +  10.616986(2*)  x  106 

(148) 

The  analysis  of  Reference  1  (Part  I  of  this  study)  demonstrated 
that  an  integration  time  of  9.6  vs  resulted  in  generally  favorable 
performance  of  the  identification  technique.  For  this  analysis, 
this  integration  is  varied  from  0.0024  vs  to  2.4  ys,  and  the  per¬ 
formance  of  the  identification  technique  is  investigated.  The 
simulation  is  set  up  to  identify  a  single  pole  system,  in  this 
case,  the  high  frequency  pole  at  s  =  -10 . 616980 (2ix )  x  106. 

The  results  of  this  simulation  are  shown  in  Table  8.  The  re¬ 
sults  indicate  that,  if  a  short  integration  time  (compared  to  the 
reciprocal  of  the  low  frequency  pole)  is  used  in  the  identifica¬ 
tion  processing,  then  the  high  frequency  pole  and  corresponding 
residue  of  H(s)  are  accurately  predicted.  The  results  indicate 
that  at  least  a  100:1  reduction  in  integration  time  from  the 
original  9.6  vs  is  required  to  effectively  isolate  the  high  fre¬ 
quency  pole  response.  These  results  suggest  that  the  integration 
time  be  less  than  l/(high  frequency  pole)  for  accurate  identifica¬ 
tion  performance. 

This  integration  time  approach  is  related  somewhat  to  the 
wide-band  processing  approach  of  Reference  2.  This  is  a  good 
technique  to  use  on  wide-band  systems  where  there  are  a  set  of 
low  frequency  poles  and  a  set  of  high  frequency  poles.  Once  the 
high  frequency  poles  and  residues  are  identified,  the  normal 
identification  procedure  is  followed  and  the  contributions  of  the 
high  frequency  poles  are  subtracted  out  from  the  system  output  be 
fore  pencil-of-functions  processing. 


In  order  to  demonstrate  how  such 
we  once  again  consider  the  system 


H(s)  = 


2.8069192  x  10l 


s  +  0 . 011550998(2* )  x  10e 


F.  SEPARATION  OF  RESPONSES 


An  important  requirement  for  the  identification  technique  is 
the  ability  to  excite  the  nonlinear  system  such  that  the  linear 
response  is  isolated  from  second-  and  higher-order  responses  and 
that  the  linear  plus  second-order  response  is  isolated  from 
third-  and  higher-order  responses.  This  requirement  impacts  both 
the  feasibility  of  implementing  the  identification  tecnnique  auu 
the  computational  complexity  involved  in  identifying  the  non¬ 
linear  impulse  responses.  This  critical  issue  is  addressed  in 
detail  in  this  section. 

The  basic  assumption  on  which  the  identification  technique 
is  founded  is  that  the  nonlinear  system  can  be  excited  in  such 
a  manner  that  the  system  response  is  linear.  Techniques  of  vali¬ 
dating  linear  operation  of  a  nonlinear  system  are  addressed  in 
detail  in  Reference  3.  Basically,  the  nonlinear  system  is  ex¬ 
cited  by  a  sinusoidal  signal  of  amplitude  A  and  a  spectral  analy¬ 
sis  of  the  system  response  is  obtained.  Amplitude  A  is  adjusted 
until  the  spectral  content  of  the  system  output  shows  that  the 
magnitude  of  second  and  higher  order  harmonic  frequencies  is 
significantly  below  that  of  the  fundamental  component.  This  pro¬ 
cedure  permits  determination  of  the  linear  impulse  response  of 
the  nonlinear  system,  hi(t).  Identification  of  hi(t)  leads,  as 
has  been  shown  previously,  to  identification  of  the  natural  fre¬ 
quencies  of  the  second  and  third-order  impulse  responses.  There¬ 
fore,  it  is  a  key  element  in  the  identification  process. 

If  the  nonlinear  system  cannot  be  excited  such  that  its  out¬ 
put  response  is  linear,  the  identification  procedure  increases 
in  complexity  but  the  second-order  impulse  response,  h2(ti,  t2) 
can  still  be  identified.  This  fact  is  demonstrated  below. 

Assume  that  the  nonlinear  system  can  be  excited  such  that 
only  y^Ct)  +  y2(t)  can  be  isolated  from  third-  and  higher-order 
system  responses.  Define  ya(t)  as 


s 


"  ~ci  ^  “  o  ji  i  j  j  B  1 . L. 

s  =  -a  t;  i  =  1 ,  .  .  .  ,L 
s  *  ;  k2  ■  1 , . . .  ,N 


(151) 


There  are  2N(N  +  1)  +  N  poles  if  L  *  N,  as  is  generally  the 
case.  Since  N  may  not  be  known,  assume  that  L  ■  1.  Then,  the 
number  of  poles  in  Ya(s)  is 


B 


2N  +  2  + 


N(N  +  1) 
2 


(152) 


These  poles  are  of  the  form 
A  i;  i  =  1,  .  .  .N 

X  i  +  A  i,  j  =  1 _ N 

-al  +  »  1  =  1 » • • »N 


-a 


1 


If  the  identification  process  is  applied  to  the  response, 
yi(t)  +  y2(t),  then  6  and  the  poles  of  Ya(s)  will  be  identified. 
By  associating  the  identified  poles  with  the  above  list,  it  will 

be  possible  to  identify  theA^,  i  =  1 . N.  The  number  of  poles 

in  the  linear  system,  N,  can  be  found  from  f} . 

This  association  will  be  done  as  follows.  In  the  list  of 
identified  poles,  there  will  be  N  pairs  of  poles  having  the  rela¬ 
tionship 


“j  *  2a£ 


a.  -  2am 
k  m 


(153) 


where  a  j  ,  a^,  a.,  ajj,  are  identified  poles.  Furthermore  there 
will  be  identified  pole  pairs  of  the  form 


an  "  "al  +  ap  (154) 

where  a-j  is  known  from  the  input.  It  is  noted  that  sufficient 

data  is  available  to  identify  the  poles  \±,  i  »  1 . N  of  the 

linear  system.  The  residues  of  the  poles  of  the  linear  transfer 
function  can  be  found  from  the  residues  of  the  poles  of  Ya(s) 
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corresponding  to  s  =  -dj..  Then  the  A^j^g  quantities  of  h2(ti,t2> 
can  be  found  in  the  normal  way. 


This  is  shown  here. 

Given  a  nonlinear  system  with  a  linear  transfer  function 
given  by 

hxc->  -  ?  fi-4-  <155> 

i-1  s  “i 
and  an  input  X(s)  given  by 
L 


X(s )  =  l 


c± 


i-1  S  +  ai 


(156) 


the  linear  response  is  described  by 

Yi 


ri<s>  ”  a  a(^)  -  rbi) 

The  second-order  response  is  given  by 

M  N  L  L  «, 

Y2(s)  *  I  Z  £  £  Ak  k  C/ 

k2-l  k^l  1*1  j*l  k1k2  1 

“i  +“j  + 


(157) 


<aJ  +  ak^  *  “kjX0!  +  aJ  +  akj  *  ak2> 

_i _ r _ i _ 

1  *  V  [(“J  +  V(“i  +  ak2) 


s  - 

( av  +  ak  ) 
1  *2 

1 

s  + 

(aj  -  V 

1 

s  + 

<°i  -  ak2> 

1 

(aJ 

■ 

+  ak2}  (ai 
\ 

i 

«i  +  aj  +  2akr 


ak  ') 
*2 


8  +  °i  +  aj 


(158) 
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If  the  response  obtained  is  y^(t)  +  y2<t),  then  Y^( s )  +  Yo(s) 
contains  poles  at  s  ■  ak2«  ^2  “  1  .,-N  whose  residues  are  of  the 

form 


The  pole  of  Y^s)  +  Y2(s)  at  s  ■  -c^  has  a  residue  given  by 


N  R  jC  j 
Z  — <3— ^ - 


J-l  “i  +  aj 


If  L  «*  N,  then  there  are  N  equations  of  the  form 


1  -  +  - - ?-4-  +  .  .  .  +  - — r~-  =  6, 


°1  +  al  ai  +  a2 


a  t  +  a 


1T  *N 


RlS  +  R2CN  + 


aN  +  al  aN  a2 


RN^N 

+  - — —  =  Bx, 

aN  +  aN  N 


or  in  matrix  form 


C1 

C1 

C1 

al + 

al 

al  + 

a_  ’  ' 

2 

*1  + 

aN 

C2 

C2 

C2 

a„  + 

a. 

a~  + 

a.  ’  '  ' 

an  + 

a7! 

2 

1 

2 

2 

2 

N 

CN 

CN 

CN 

aN  + 

~*1 

aN  + 

a2 

aN  + 

aN 

N 


N 


or  W  W-  [»]  in  matrix  form. 

For  linear  independence,  it  is  necessary  that 
det  f  0 
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(159) 


(160) 


(161) 


vs*fTw »47*r^ i" r  A* 


1 1 


This  has  been  shown  in  Reference  3  provided 
f  aj  for  any  i,  j  -  1,...,N;  i  f  j 

ai  1  Q j  for  any  i,  j  »  1, . . . ,N;  i  ^  j 

as  is  the  case  for  the  identification  technique. 

Although  the  procedure  is  more  complicated,  the  linear  and 
second-order  impulse  responses  can  be  identified  even  if  the 
second-order  response  cannot  be  isolated  from  the  linear  response. 
The  practical  application  of  this  procedure  may  be  complicated  by 
the  need  to  determine  6.  This  is  a  potential  numerical  accuracy 
problem.  However,  in  theory  at  least,  the  linear  response  need 
not  be  isolated  from  the  second-order  response. 

Another  complication  in  the  identification  procedure  arises 
if  the  second-order  response  cannot  be  isolated  from  the  third- 
order  response.  Assume  that  the  linear  impulse  response  of  a 
nonlinear  system  has  been  identified.  If  the  system  cannot  be 
excited  such  that  the  response  yi(t)  +  y2(t)  is  obtained,  then 
the  identification  technique  will  use  the  response  yj.(t)  +  y2(t) 

+  y3(t) ,  where  it  is  assumed  that  fourth  and  higher  order  re¬ 
sponses  are  negligible  compared  to  third  order.  For  convenience, 
we  define 

yb(t)  =  y2(t.)  +  y3(t)  (162) 

and 

Yb(8)  =  Y2(s)  +  Y3(s)  (163) 

where  y-^t)  has  been  subtracted  out  using  knowledge  of  Hi(s)  and 
hi(t).  The  poles  of  Yt>(s)  are  given  by 


s  = 

ak  +  ak  +  ak  ; 
K1  k2  *3 

;  k1  =  1 , . . . 

,J;  k2  =  1,. 

.  M ;  k1  =  1,  .  . 

.  .  ,N 

s  » 

-  ai  *  ak  + 

•3 

;  i  ,  k  =  1 , .  . 

•t  L  t  ^  C  ^  »  * 

.  . ,  N 

s  = 

'  ai  "  aJ  ‘  “k; 

i,j,k  -  1,. 

.  •  t  L 

s  *=  a^  +  a^  ;  =  1,...M;  kg  *  1,...,N 


s  «  -  ai  +  a^  ;  i  =  1 , .  .  .,  L;  kg  = 
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s  “  ~  ai  "  V  i,J  "  1>L 

s  =  -  ai  +  ak  +  ak  ;  i.kg  -  1,N;  k2  =  (164) 


All  of  these  quantities  are  known,  since  it  was  assumed  that 
h^(t)  has  been  identified.  The  problem  is  to  identify  the  Afcik2 
of  ho(ti.to).  The  poles  of  Yk(s)  which  contain  information  about 
the  are  Siven  by 

s  =  a.  +  a.  ;  k1  *  k„  * 

K1  k2  1  * 

s  =  -  ai  +  ;  i  =  1,  .  .  .  ,L;  k2  *  1 ,  .  .  .  ,N 

2 

s  a  '  ai  "  aj >  i,J  =  1  * • • •  •  L 

The  poles  of  Y2(s)  that  correspond  to  poles  of  Y3(s)are  those 
for  which 


ak  +  ak  +  ak  =  ak  +  ak 
K1  k2  k3  k4  5 

k^  -  1 ,  .  .  .  ,  J ;  k2  =  1 , . . . , M ;  k^  *  k^  *  1 , .  .  .  , M 

k5  *  1 , . . . ,N 


and 


+  a, 


k2,k4 


N,  k . 


1, 


M 


(166) 


These  correspond  to  poles  of  the  form 
Xj  J  " 

XA  +  Xj  £,J  -  1, . . .  ,N 

-  “i  +  Xj  i.J  =  1, - N  (167) 

If  the  residues  of  these  poles  are  known,  then  a  set  of  equ¬ 
ations  exists  that  involves  pairs  of  A^ko  and  C,  .  quantities. 

J-  *  *1*2*3 

The  key  question  is  whether  the  A^k2  and  ckikok3  quantities 
can  be  determined. 
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Suppose  we  consider  the  pole  at  s  =  X^  +  X2-  The  residues 
involve  the  following  quantities 


A12’  C1 , 2 , 1 ‘  C1 , 2 , 2 ' 


'1.2.N 


This  implies  that  there  are  N  +  1  unknown  coefficients  in  the 
residue  equation  for  the  pole  at  s  ■  X^  +  Xg.  It  is  necessary 
to  determine  if  these  quantities  can  be  solved  for  using  the 
residues  of  Y2(s)  +  Y3(s).  It  is  shown  in  Reference  3  that  the 
ck3kok3  quantities  of  the  form  Cjj^,  i,j  *» 

i  f  j,  are  identified  from  the  poles  of  Y3(s)  given  by  s  *  X^ 

+  X j  +  X.  ;  i.j.k  =  1 . N;  i  f  j ,  j  ^  k.  Since  these  poles  are 

unique  to  Y3(s),  the  above  j ji  ,  and  quantities  can 

be  determined  in  the  usual  manner.  J0nce  these  are  known,  the 
Ajj  quantities  are  found  directly  from  the  residue  at  s  =  Xi  +  X2. 
Similarly,  all  the  Akik2»  ki,  k2  *  ki  f  k2 ,  can  be  found 

from  the  poles  of  Y2<s)  +  Y3(s)  at  s  =  Xki  +  Xk2* 

In  Reference  3,  it  is  shown  that  the  Ckik2k3  are  identified 
from  residues  of  the  poles  at 


s  = 

3xi 

i  =  1, . . . ,N 

s  = 

2Xi  ♦  Xj 

i , j  *  i  f  j 

s  - 

Xi  +  Xj  +  Xk 

i.J.k  »  1.....N;  if 

k 

s  = 

'  ai  +  Xj  +  Xk 

i  fixed,  j  f  k;  J,k  =  1,, 

.  .  .  ,N 

s  = 

"  “i  +  Xj 

j  -  1, . .  .  ,N 

(168) 

All  of  these  poles  can  be  used  to  identify  the  Ckikoks  as  is 
normally  done  except  for  the  poles  at  s  s  -  a,  +  XV  JWe  now 
consider  the  portion  of  the  response  Y2(s)  +  Y3(s)Jdue  to  the 
pole  at  s  =  ~ai  +  Xj.  The  unknown  Ckik2k3  quantities  are  of  the 
form  Ck^mn  where  ki>N,  m>N ,  n<N  and  m  and^n  are  such  that 


are  such  that 


ln  =  Xl- 

ni  1 


(169) 


The  unknown  A 


k  k  quantities  are  of  the  form  A^ ,  g, 

12  1  1 

and  a.  +  a.  =  X«.  There  are  N  pairs  of  (a 


where  k^  >  N  and  ak  +  aj  =  X£.  There  a 

such  that  am^  +  aR^  =  X^  The  portion  of  the  third  order  response 
due  to  the  pole  s  *  -a^  +  X^  is  given  by 


N 

kl  N 

V 
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■'UH? 


i  cWl  j  i  (X  A 

Y3c(=)*ki.U  i-h-i[  \rvx^V%?j 

_  /  aJ-*-a.+2am  \ 

/  \  Ck«in«nrt  L  L  /  1  J  2  \ 

/ _ 1 _ \  +  -  }  2-~  I  £  -3  l(a . +a  )(a,+a“  )J 

l  ai+ctj+an>1+an1y  J  *1  %  i-1  j  =  1|_  \  1  m2  2/ 


!  \“|  CklmNnM  |  j 

vvv\)  *  ‘  ai+‘ki  1=1  J=1 


°i+V2a*M 


[-3  ^v^vv/y^vx n» 


S+Oj-X 


(170) 


The  portion  of  the 
given  by 


second-order  response  it  s  =  -  «i  *  'i  ls 


Ws)  •  * 


L  ,  /  AklN  \  ^  1 

*ml  (d1  7  a^T  +  ■  ’  •  +  (“i  +  aN 


(171) 


o  ,  „  p  Quantities  and  n2  unknown  Aklk2 

There  are  2n3  unknown  cklk?^  t  Qf  the  response, 

quantities  in  the  Y2c(s)  «■  ** 

+  =  _  a,  +  a.  is  defined  to  be  Q1.  The 

functional8 form  of  ^residue  Ian  hi  written  as 
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c  c 

J  k-m.n- 

3  l  — -  - ■  ft  +  22  2  ft  +. 

kj-N+1  al  +  akl  1  a2  +  akx  2 

M  Akxl  Akj2 

kj-N+1  al  +  ak;L  Yl  a2  +  ak2  Y2 

This  can  be  reduced  to 


klmNPy 

aN  +  *k1  N 


k.N 

.  +  - = -  v„ 

aM  +  a.  'N 
N  kj 


(172) 


kl=N+l 


( !SvvLi  +  VYi\ /VrtV iVji + 

V  /\  a»*\  / 

/3Ck,m„n  8  ,  ,  +  Ak<NY,\]  3  P^k-m,  n  8  3Ck,ro_n  8 


+  f  klmNnN  N  klN  N\  ,  E  kimiPl  1  +  kim2n2  2 

\  aN  +  %  J  kl=M+1  Q1  +  akl  “2  +  akl 


+  +_0_0_N=  @1 

J  kn 


(173) 


Let  Ej  =  then  rewrite  equation  (173)  as 


'  °i  *\~) 1  v  °»  +  n  ] 


J  3Cklm  n  3Cklm  n  8 

l  - ±_i_i  g  +...+  - *-  a  e, 

k^=M+l  ®1  ak;L  1  aN  akx  1 

Furthermore,  equation  (173)  can  be  rewritten  as 


(174) 


2N  F 

^  a  = 

i-1  a,  +  f. 


(175) 
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where  F.  »  C,  B  +  A  r 

j  kimini  i  kiici 


for  ki  which  have  nonzero 

Ck1tnlni  and  V°  for  1>N 


A  similar  equation  to  (175)  is  obtained  for  each  value  of 
ai >  If  L«2N2,  then  the  identification  technique 

generates  the  equations 


2N 

Z 


2IT  F. 

I  - J — : 

j-1  a1  * 

2 


J  +  1  a2 


J. 


+  f 


a  “  62 


(176) 


2N* 

Z  • 
j=i  a2N: 


2  +  f 


a  "  92N2 
j 


The  set  of  equations  in  (176)  was  previously  shown  to  be  linearly 
independent.  Therefore,  the  set  of  equations  provides  a  unique 
solution  for  the  Fj  quantities. 

There  still  remain  N  pairs  of  and  Akik2  terms  which 

have  not  been  identified.  These  are  of  the  form 

3CV  m  n  +  K  lci 


These  quantities  must  be  separated  to  completely  identify 
h2(t,,t2)  and  h^tj.  tg.tg).  There  is  a  need  to  find  a  way  to 
ir 


sepafatfily  identify ’ t^ese  A 


klk2 


and  C 


klk2k3 


quantities. 


quan- 


Mike  poles 


The  only  source  for  unique  identification  of  the  Aj 
titles  is  the  residues  of  the  poles  at  s  *  -a,  -  a., 
are  unique  to  Y£(s)  and  the  residues  involve  only  the  unknown 
Akik2  quantities.  The  problem  is  to  demonstrate  whether  or  not 
these  residues  can  be  used  to  generate  a  set  of  linearly  inde¬ 
pendent  equations  to  permit  solution  for  the  Ajtlj59.  Consider  the 
pole  at  s  *  -2ar  - --  ' -  1  * 


The  residue  is  of  the  form 
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i  j  ?  »  r  i  1 
t°rTTi)  k^-i  k2-i  k*ik2  (ai +  ak2} 


i  ;  ?  i  ■ 

(20tl  +  X2}  k-'^l  k2-l  Ak"lk2  (al  +  ’ak2) 


*  (2«!  *  X2,  N,^  *  V*2  (“1  *  %) 


where  k.',  k-",  ...  k.  correspond  to  those  values  of  k.  for 
which  are  unknown.  This  equation  involves  n2  unknown  Akik2 

quantities?  Analysis  has  not  been  able  to  show  that  this  set 
of  equations  is  linearly  independent  or  dependent.  The  equations 
of  interest  for  linear  independence  become 

1  M  Ci  +  !  2N  C1 

ial  +  al  i=l  “l  +  ai  2ctl  +  a2  i=N+l  “l  +  ai 


• • •  +  2a  +  aN  1  2  a  +  at  =  0 
1  w  i-n  -N+l  1  X 

N  c  2N  C 

_ x  r  _ ± +  _ ± _  £  _ = — 

20g  +  a^  ias^  a2  +  a^  2o2  +  a2  i=N  a2  +  a^ 


1  w  vi 

+  1  r  _ m  Q 

2a„  +  aN  9  a,  +  a. 
i  N  i-N^-N+l 


JctN  +  al  i-1  aN  +  aN  2ct2  +  a2  i-N+1  aN  +  al 


+  +  x  y  - ± —  s  o 

2aN  +  aN  2  °N  +  ai 

N  N  i»N-N+l 


(177) 
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If  this  set  of  equations  implies  C.*0,  i“l,  ....  N  ,  then  the  set 
of  equations  involving  the  A^j^  i®  linearly  independent.  Succes¬ 
sive  solution  of  these  simultaneous  equations  does  not  produce  a 
factorable  polynomial  to  demonstrate  independence.  The  resultant 
equation  cannot  be  solved  and  linear  dependence  or  independence 
cannot  be  shown.  A  similar  situation  exists  for  the  unknown 
Ckik2k3  quantities  and  the  poles  at  s  «  -ai  -  -  <*k  that  are 

unique  to  Y3(s). 

This  analysis  has  failed  to  demonstrate  that  there  is  no  need 
to  isolate  the  second-order  system  response  from  the  third-order 
response.  Therefore,  it  is  concluded  that,  for  a  practical  im¬ 
plementation  of  the  identification  technique,  it  is  necessary  to 
excite  the  nonlinear  system  so  that  the  third  order  and  higher 
order  responses  are  negligible  compared  to  the  linear  and  second 
order  response. 

G.  ALTERNATIVE  IDENTIFICATION  PROCESSING  ALGORITHMS 

The  identification  technique  described  in  this  report  is  based 
on  the  pencil  of  functions  approach  to  linear  system  identifica¬ 
tion.  The  first  step  in  the  nonlinear  identification  process  is 
the  identification  of  the  poles  and  residues  of  the  linear  system 
transfer  function.  The  pencil-of-functions  approach  is  well- 
suited  to  accomplishing  this  identification.  The  second  step  in 
the  nonlinear  system  identification  process  is  the  identification 
of  the  residues  of  the  poles  of  Y2(s).  These  poles  of  Y2(s)  are 
known  once  the  linear  transfer  function  is  identified.  The 
pencil-cf-functions  approach  is  still  used,  but  it  must  be  noted 
that  other  potential  methods  of  identifying  the  residues  of  the 
poles  of  Y2(s)  exist  that  might  alleviate  some  of  the  computa¬ 
tional  complexity  associated  with  the  pencil-of-functions  approach. 

The  basic  problem  at  this  point  is  to  identify  the  R^  quanti¬ 
ties  in  the  equation 


y2(t)  = 


,-ckt 


(178) 


where  y2(t)  is  the  second  order  response  of  a  nonlinear  system. 

In  this  equation,  the  are  known,  since  they  are  related  to  the 
poles  of  the  linear  transfer  function  of  the  system.  A  sampled 
time  history  of  y2(t)  is  obtained  via  measurement  of  the  output  of 
the  nonlinear  system.  The  objective  then  is  to  use  this  infor¬ 
mation  to  evaluate  the  R^ . 


-  ■  j 


A  detailed  review  of  several  candidate  algorithms  for  solving 
this  problem  has  been  addressed  in  Reference  6.  These  identifi¬ 
cation  approaches  include  the  least-squares  method,  orthonormal 
least  squares  method,  equality  of  derivatives  method,  equality  of 
integrals  method  and  the  generalized  integrated  squared  error. 
Details  of  these  approaches  are  given  in  Reference  6  and  are  not 
repeated  here. 

The  basic  approach  used  in  many  of  these  methods  is  to  sample 
the  time  function  (in  this  case,  ygtt))  M  times  where  M  >  8  [6  is 
the  number  of  poles  in  Y2(s)  or  natural  frequencies  in  y2(t)}  • 
Then,  each  technique  attempts  to  minimize  an  error  function  to 
determine  an  "optimum"  set  of  Rfc  coefficients.  In  general  these 
techniques  require  inversion  of  a  B  x  B  matrix  to  determine  the 
Rk  coefficients.  The  matrix  entries  involve  functions  of  the 
e-skt  and  in  this  sense  are  very  similar  to  the  technique  used  in 
the  pencil-of-functions  approach.  The  advantage  of  the  pencil-of- 
functions  method  is  that  no  approximations  are  used  as  is  the 
case  with  these  overdetermined  system  approaches  (M  >  B).  Because 
of  this  advantage  and  the  requirement  that  a  B  x  B  matrix  be  in¬ 
verted  by  these  other  identification  techniques,  the  pencil-of- 
functions  approach  appears  to  be  as  good  a  candidate  for  this 
nonlinear  system  identification  technique  as  the  others  described 
in  Reference  6. 
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